Skip to main content

Phthalates impact on the epigenetic factors contributed specifically by the father at fertilization



Preconception exposure to phthalates such as the anti-androgenic dibutyl-phthalate (DBP) impacts both male and female reproduction, yet how this occurs largely remains unknown. Previously we defined a series of RNAs expressly provided by sperm at fertilization and separately, and in parallel, those that responded to high DBP exposure. Utilizing both populations of RNAs, we now begin to unravel the impact of high-DBP exposure on those RNAs specifically delivered by the father.


Enrichment of RNAs altered by DBP exposure within the Molecular Signature Database highlighted cellular stress, cell cycle, apoptosis, DNA damage response, and gene regulation pathways. Overlap within each of these five pathways identified those RNAs that were specifically (≥ fivefold enriched) or primarily (≥ twofold enriched) provided as part of the paternal contribution compared to the oocyte at fertilization. Key RNAs consistently altered by DBP, including CAMTA2 and PSME4, were delivered by sperm reflective of these pathways. The majority (64/103) of overlapping enriched gene sets were related to gene regulation. Many of these RNAs (45 RNAs) corresponded to key interconnected CRREWs (Chromatin remodeler cofactors, RNA interactors, Readers, Erasers, and Writers). Modeling suggests that CUL2, PHF10, and SMARCC1 may coordinate and mechanistically modulate the phthalate response.


Mediated through a CRREW regulatory network, the cell responded to exposure presenting stressed-induced changes in the cell cycle—DNA damage—apoptosis. Interestingly, the majority of these DBP-responsive epigenetic mediators’ direct acetylation or deacetylation, impacting the sperm's cargo delivered at fertilization and that of the embryo.


It is now well-known that sperm delivers an entire set of extra-chromosomal components, including RNAs, at fertilization ([1], reviewed in [2]). We defined a series of paternal-provided RNA elements (REs, exon-sized sequences) that are enriched at least fivefold above the oocyte and delivered at fertilization, providing a unique set of RE-containing RNAs (RE-RNAs) while markedly enhancing those present in the oocyte [3]. To date, several of these sperm RE-RNAs have been shown to respond to exposures reflective of lifestyle [4,5,6,7]. Some have now been implicated in offspring phenotype [6, 8,9,10,11], highlighting the importance of understanding their role in development.

Phthalates are endocrine disruptors widely used in consumer products [5, 12,13,14], such as vinyl plastics, personal care products, and some medication coatings [15,16,17,18]. To date, it is known that phthalate exposure in males has an adverse impact on semen and embryo quality, as well as time to pregnancy [19,20,21,22], but the mechanism(s) remains obscure. To begin to address this gap, we defined a series of sperm REs that respond to phthalates [5] using the dibutyl-phthalate (DBP) Inflammatory Bowel Disease (IBD) mesalamine crossover cross-back model (reviewed in [17]). To evaluate the impact of high-DBP exposure, we recruited men taking one of two formulations of mesalamine; one was encapsulated in a DBP-containing coating, and one was without DBP in the coating. Sperm RE-RNAs [5] were isolated and compared to those observed in the oocyte and zygote. Two paternal provided classes were defined from REs present in the zygote. Those paternal provided REs enriched ≥ fivefold compared to the oocyte [3], and excluding those defined here as fivefold enriched [3], those paternal REs, ≥ twofold enriched compared to the oocyte. Through RE expression, we show the mechanistic impact of high-DBP exposure acting through epigenetic modifiers and how they affect those RE-RNAs paternally provided to the oocyte at fertilization.


Sperm RNAs are known to respond to environmental exposures [4,5,6]-like dibutyl-phthalate (DBP), an endocrine disruptor found in some medications, including the coating of Asacol, whose active ingredient is mesalamine used to treat Inflammatory Bowel Disease (IBD) [5, 15, 22]. At the recommended maximal Asacol dosage, DBP exposure from the coating exceeds the Environmental Protection Agency (EPA) reference dose for a 150-pound individual by 300–700% (reviewed in [5, 23]) based on the DBP primary urinary metabolite, monobutyl phthalate (MBP). Men on Asacol had MBP urinary concentrations 1,000 times higher than the median male concentration reported in the National Health and Nutrition Examination Survey (NHANES [24] within the general United States population [23]. However, studies also indicate that lower level environmental background exposures to DBP from personal care and consumer products may impact semen and embryo quality, and time to pregnancy [19,20,21,22]. Analysis of those sperm RNA Elements (REs, exon-sized RNA fragments) altered in response to high-DBP exposure from using DBP-coated mesalamine, Asacol, compared to the non-DBP coated mesalamine, Pentasa, has begun to define DBP exposome pathways that impact sperm RE-containing RNAs (RE-RNAs) [5].

As summarized in Fig. 1A, REs responsive to high-DBP exposure [5] modifying the male contribution at fertilization was considered. Men who were on non-DBP coated mesalamine (e.g., Pentasa) transitioned to high-DBP coated mesalamine (Asacol) (referred to as baseline to crossover; B1H) in the baseline non-DBP (B1HB2) study arm as well as men transitioning from non-DBP coated mesalamine (e.g., Pentasa) back to high-DBP coated mesalamine (Asacol) (referred to as crossover to crossback; BH2) in the high-DBP (H1BH2) study arm, were considered. Comparison of REs responding to DBP withdrawal was from the men starting on high-DBP coated mesalamine transitioning to non-DBP coated mesalamine (baseline to crossover; H1B) and men transitioning from high-DBP-coated mesalamine back to non-DBP-coated mesalamine (crossover to crossback; HB2).

Fig. 1
figure 1

Study analysis design. A Briefly, dibutyl-phthalate (DBP) responsive RNA Element (RE)-containing RNAs (RE-RNAs) from each crossover–crossback segment and paternally provided set were identified. B Ontology was assessed by the Molecular Signature Database (MSigDB), and shared DBP responsive and paternally provided RE-RNAs determined. C Shared RE-RNAs were evaluated for presence in overlapping DBP responsive and paternally provided MSigDB enriched gene sets. CRREWs (chromatin remodeler cofactors, RNA interactors, Readers, Erasers, and Writers) RE-RNAs were identified, and their presence within overlapping DBP responsive and paternally provided MSigDB gene sets was assessed. D Gene network to visualize CRREW biological process interactions was generated. H1B; high-DBP (baseline visit) to background DBP (crossover visit), BH2; background DBP (crossover visit) to high-DBP (crossback visit), B1H; background DBP (baseline visit) to high-DBP (crossover visit), HB2; high-DBP (crossover visit) to background DBP (crossback visit)

Biological pathways altered by DBP exposure

The REs defined from this series of high-DBP exposures or withdrawals were associated with biological pathways within the Molecular Signature Database (MSigDB) using the DBP responsive RE-RNAs summarized in Fig. 1B (Additional file 2: Table S1) as input. These included all positively or negatively correlated RE-RNAs, or the entire set of correlated RE-RNAs described as altered between each study visit, e.g., B1H (high-DBP exposure). For each set of genes, enrichment (Fig. 1B) identified five major biological processes; cellular stress, cell cycle, DNA damage response, apoptosis, and gene regulation (Fig. 2A, Additional file 3: Table S2). Each enriched biological process included specific MSigDB gene sets within both study arms (B1HB2 and H1BH2).

Fig. 2
figure 2

Biological pathways enriched in response to dibutyl-phthalate (DBP). A DBP responsive enriched biological pathways were identified using the Molecular Signature Database (MSigDB). The total number of unique RE-RNAs and enriched gene sets is provided. Enriched gene sets were either unique to DBP study arm (B1HB2 or H1BH2) or within each study arm. B1HB2 provided the background DBP study arm baseline to crossover to crossback, H1BH2 provided the high-DBP study arm baseline to crossover to crossback. B DBP-responsive RNA Elements (REs) within the enriched biological pathways that share paternally provided RE-containing RNAs (RE-RNA) were identified. Paternally provided: REs determined as fivefold paternally enriched or twofold paternally enriched. Node color indicates the following; Orange: enriched biological processes, Yellow: RE abundance increases following DBP exposure, Pink: RE abundance decreases upon DBP exposure, Green: RE abundance increases upon DBP withdrawal, Blue: RE abundance decreases upon DBP withdrawal. Edge color indicates the total number of DBP responsive REs (1–40 REs) with an overlapping paternally provided RE-associated gene name moving between nodes. Scaling for edge color is continuous with 1 RE attributed to a Dark Blue color and 40 REs attributed to a Red color. Edge line type indicates the following; solid: background DBP (baseline visit) to high-DBP (crossover visit) (B1H), parallel lines: high-DBP (crossover visit) to background DBP (crossback visit) (HB2), dash: high-DBP (baseline visit) to background DBP (crossover visit) (H1B), dots: background DBP (crossover visit) to high-DBP (crossback visit) (BH2)

The majority of enriched gene sets were associated with gene regulation (191/392 gene sets), indicating a large proportion of DBP-responsive RE-RNAs function as either transcription factors (TFs) or CRREWs (Chromatin remodeler cofactors, RNA interactors, Readers, Erasers, and Writers [4]). DBP-responsive TF binding sites were identified, and the number of unique TFs assigned to these binding sites is summarized in Additional file 4: Table S3A. However, there was no statistical significance in the number DBP responsive TF encoded RE-RNAs (Additional file 4: Table S3A), which suggested other modulators of gene expression may dominate, and CRREWs were considered (Fig. 1B and D, Additional file 5: Table S4). The proportion of CRREWs was higher than expected (Additional file 4: Table S3B). Accordingly, their potential as modulators among the various enriched biological processes was examined (Fig. 2A). In total, 119 CRREWs were specific to the analysis within the B1HB2 comparisons, while 60 CRREWs were specific to those within the H1BH2 comparisons. Fifty CRREWs were represented in the analysis of both study arms (Additional file 5: Table S4). Together this suggests that a start arm's initial drug coating (i.e., high-DBP coated mesalamine or non-DBP coated mesalamine) impacts these biological processes through a unique set of CRREW modulators. In total, 50 DBP-responsive CRREWs were shared between both study arms modulating these biological processes.

DBP-responsive and paternally provided RE-RNAs

Consideration was given to whether the DBP-responsive RE-RNAs and enriched pathways were present, enriched, and/or unique in those provided by the father at fertilization. Enriched vs unique paternal REs were defined by the presence, or lack of the RE in the oocyte, with an RPKM < 2 considered absent, as this is the lower threshold in which RE presence exceeds experimental error. Two types of paternally provided REs, and their associated gene names (RE-RNAs), were examined. The first comprised the 289 REs (from 206 RE-RNAs) previously defined as paternally enriched (fivefold paternally enriched, Fig. 1A) [3] from a non-IBD population exposed only to background levels of DBP. The second included an additional set of 250 REs (from 93 RE-RNAs), not including those previously defined, that appear to have at least a twofold enrichment in sperm compared to the oocyte (twofold paternally enriched, Fig. 1A, Additional file 6: Table S5). Both paternally provided sets were identified from those observed at a level of at least 10 RPKM in the zygote. A total of 83 REs (5 × paternally enriched: 75 REs [3], 2 × paternally enriched: 8 REs (Additional file 6: Table S5)) were specific to the sperm. In addition, as previously stated, an RPKM > 10 in the zygote was a requirement for the paternally provided REs. This leads to the possibility of the RPKM abundance in the zygote exceeding the combined RPKM in the sperm and oocyte. Only 2% of the paternally provided RE-RNAs (11/539 RE-RNAs) > 5 RPKM were more abundant in the zygote than in sperm. Comparing these paternally provided RE-RNAs to those defined as DBP responsive defined a series of RE-RNAs that were both DBP responsive and paternally provided (Figs. 1B, 2B).

Using Cytoscape [25] to aid visualization, the 132 paternally provided DBP-responsive RE-RNAs identified in Fig. 2B represented a significantly larger overlap between data sets than expected (p < 4.358e-35, representation factor = 3.0). Of these RE-RNAs, a significant proportion was CRREWs (45 CRREWs; Fig. 1C, Table 1, Additional file 7: Table S6). From this list, 86 RE-RNAs, including 38 CRREWs (Table 1), were within the DBP-enriched biological processes, including GOMF chromatin binding and GOBP cellular response to stress (Figs. 1B, 2B). Most of the DBP responsive CRREWs were specific to the B1HB2 or the H1BH2 comparisons. The chromatin remodeler cofactors PHF10, CUL2, ATXN7, PHC3, SMARCC1, and EYA3, were identified (Table 1), indicating these CRREWs likely modulate the DBP response between these processes. Of these, CUL2, PHF10, and SMARCC1 were visually full-length in all DBP responsive and paternally provided samples that passed the Transcript Integrity Index (TII) for identification of samples of similar RNA quality [26] (Table 1, Additional file 1: Figure S1). Their representation in the shared enriched gene sets within each biological process supports the essential role of these three CRREWs in the DBP response.

Table 1 Dibutyl-phthalate (DBP) responsive Chromatin remodeler cofactors, RNA interactors, Readers, Erasers, and Writers (CRREWs)

As summarized in Table 2, the paternally enriched DBP-responsive RE-RNAs, including ACSM3, PSME4, CCDC7, NUP98 and CAMTA2, responded similarly throughout the study, increasing or decreasing in abundance when exposed to or withdrawn from high-DBP. CAMTA2 and PSME4 were full-length (Additional file 8: Table S7), suggesting post-fertilization activity. The corresponding MSigDB gene sets for the aforementioned DBP responsive and paternally provided CRREWs, along with the non-CRREW RE-RNAs CAMTA2 and PSME4, are highlighted in Fig. 3 and Table 3.

Table 2 Dibutyl-phthalate (DBP) responsive and paternally provided RNA Element (RE)-containing RNAs (RE-RNAs)
Fig. 3
figure 3

Interconnected gene network of paternally provided dibutyl-phthalate (DBP) responsive RNA Elements (RE)s. The gene network highlights those relationships corresponding to the full-length RE-containing RNAs (RE-RNAs) CAMTA2, EFCAB6, PSME4, SMARCC1, PHF10, RARA, FBRS and CUL2. Paternally provided RE-RNA is either fivefold or twofold paternally enriched. Key genes are highlighted based on node color. Dark blue node borders indicate fivefold paternal enrichment, while a green border indicates twofold paternal enrichment. Node color and shape are as follows; pink squares: CRREW (Chromatin remodeler cofactor, RNA interactor, Reader, Eraser, and Writer) RNA, blue square: key RNA that is not a CRREW, orange circle: major enriched biological pathways, bright yellow circle: indicate functions related to acetylation/deacetylation, methylation/demethylation and ubiquination/deubiquination: light yellow circle; specific process related to acetylation/deacetylation, methylation/demethylation and ubiquination/deubiquination, light pink circle: protein complexes that include at least one key gene, purple circle; enriched MSigDB geneset, grey circle; gene function not assigned by an MSigDB geneset, blue circle; MGI phenotype

Table 3  Enriched biological processes of key dibutyl-phthalate (DBP) responsive, paternally provided RNA Element (RE)-containing RNAs (RE-RNAs)

A total of 46 of the 132 DBP responsive and paternally provided RE-RNAs were not within any of the enriched MSigDB gene sets related to the biological processes highlighted in Fig. 2. These 46 RE-RNAs were evaluated separately using EnrichR to identify gene function (Additional file 9: Table S8) and determine whether they reported functions within these biological processes. This highlighted multiple full-length CRREWs. First is the histone deacetylase complex (HDAC) [27, 28] regulated EFCAB6, a chromatin remodeler cofactor associated with androgen receptor (AR) signaling, proteolysis, and transcriptional regulation. It further highlighted the chromatin remodeler cofactor RARA, which functions in histone methylation [29, 30], is involved in apoptotic cell clearance, the positive regulation of the cell cycle, cellular response to estrogen and hormone stimulus, and the negative regulation of transcription by RNA polymerase II (Additional file 9: Table S8). Interestingly, the chromatin remodeler cofactor FBRS had no reported ontologies within EnrichR. However, it is part of the RING2-FBRS replication-independent complex involved in histone modification [29, 31, 32].


Utilizing the DBP-responsive sperm RE-RNAs [5], fivefold paternally enriched RE-RNAs [3], and twofold paternally enriched RE-RNAs, we have begun to frame the effect of exposures on what is paternally provided at fertilization. Significantly more than expected RE-RNAs corresponding to CRREWs (38 DBP responsive, paternally provided CRREWs within shared enriched MSigDB gene sets) were identified. This included three full-length RE-RNAs, consistent with the view that they act as mediators between pathways in response to exposure. This increase in the proportion of CRREWs but not TFs mimics what we previously reported within a series of body mass index (BMI) responsive RE-RNAs [4]. Consistent with the above, if they play a role during protamine replacement and syngamy as the blastocyst begins to form, one would expect a more significant proportion of CRREW RE-RNAs than TFs remaining in the mature sperm.

DBP impacts 86 paternally provided RE-RNAs (Figs. 2B, 3) ontologically linked to cellular stress, cell cycle, apoptosis, DNA damage response, and gene regulation. While most of these RE-RNAs appear responsive to DBP in either the H1BH2 or B1HB2 exposure comparisons, a subset is responsive irrespective of exposure duration or time removed from high-DBP. Three full-length CRREWs (CUL2, PHF10, and SMARRC1) and two full-length fivefold paternally enriched non-CRREW RE-RNAs (CAMTA2 and PSME4) were identified. These five RE-RNAs respond to the addition or removal of high-DBP irrespective of the study arm and are within a series of biological processes, as represented in Fig. 3.

This analysis identified a highly complex, interconnected gene network (Fig. 1D) reflecting DBP responsive and paternally provided RE-RNAs. To more readily interpret this gene network, the focus was given to the full-length DBP responsive paternally provided RE-RNAs highlighted in Fig. 3. Here, SMARCC1 and PHF10 (top of Fig. 3) function as part of the SWI/SNF complex (middle left of Fig. 3) to enhance the transactivation of AR [33,34,35] in direct opposition to the AR transcriptional repression by the chromatin remodeling cofactor EFCAB6 (Additional file 9: Table S8, bottom of Fig. 3) [27, 28]. Interestingly, DBP has anti-androgenic activity (reviewed in [13]); however, AR was not found as directly responsive to high-DBP exposure or withdrawal [5], which is consistent with the lack of an in vitro interaction between DBP and AR [36]. In response to DNA damage, the SMARCC1 and PHF10 containing SWI/SNF complex (middle left of Fig. 3) are known to accumulate, likely enabling transient chromatin accessibility to DNA-binding and DNA damage response proteins [37]. As expected, PHF10 and SMARCC1 were within a number of the same enriched gene sets (Table 3, Fig. 3), including GOMF histone binding alongside PSME4, and GOMF transcription regulator activity and GOBP positive regulation of transcription by RNA polymerase II with CAMTA2, indicating shared functions. CAMTA2, while not a CRREW, is a transcriptional activator that associates with class II HDACs to negatively modulate topological associated domains (TADs) [38, 39]. However, the proteasome component PSME4 acts to recognize acetylated histones, promoting histone degradation during spermatogenesis and the DNA damage response [40, 41]. The chromatin remodeler cofactor RARA regulating transcription in a ligand dependent manner (Additional file 9: Table S8). RARA functions in response to estrogen stimulus (Additional file 9: Table S8), is involved in H3K4 methylation [30] as part of a heterodimer and induces histone deacetylation when the heterodimer associates with specific multiprotein complexes [42]. These relationships begin to highlight how the DNA damage response may feedforward regulating gene expression (center of Fig. 3). Together with PSME4, CUL2 was within the enriched DNA damage response gene set Dacosta UV response via ERCC3 dn (Table 3, Fig. 3). As part of the E3 ubiquitin ligase complex (middle left of Fig. 3), CUL2 alongside BAF250, elongin C and ROC1 ubiquitinate histone H2BK120 aiding in SWI/SNF complex H3K4 trimethylation [43, 44]. CUL2 further enables the interaction of VHL with elongin B and elongin C to form the E3 ubiquitin ligase complex to recruit VHL to HP1 chromatin [45, 46]. In addition, the CUL2 containing E3 ubiquitin ligase has been identified as important in Adenovirus inactivation of a DNA damage response [47]. It is integral to the progression of G1 to S and the S-phase-dependent DNA damage response [48]. While not identified within the enriched MSigDB gene sets from Additional file 3: Table S2, RARA participates in the regulation of the cell cycle and apoptotic cell clearance (Additional file 9: Table S8, bottom of Fig. 3). This indicates a potential interaction with CUL2 within the apoptotic gene sets Graessman apoptosis by doxorubicin dn and GOBP programmed cell death (Table 3, Fig. 3). These relationships highlight the potential cooperation between these CRREW complexes and how a DBP-induced DNA damage response may impact the cell cycle, leading to its arrest and eventually, cell death.

CUL2 and PSME4 were also within REACTOME cellular responses to stimuli (Table 3, Fig. 3). While little is known about the function of FBRS, it is part of the RING2-FBRS complex (bottom right of Fig. 3), a type of Polycomb group (PcG) complex [29, 31, 32] that acts as a transcriptional activator of mesoderm differentiation, and a regulator of H2AK119ub1 levels [49]. PSME4 and SMARCC1 are within the gene sets REACTOME RNA polymerase II transcription, Senese HDAC3 targets up, and GOBP chromosome organization (Table 3, Fig. 3). With DBP-responsive RE-RNAs represented within each enriched biological process, an interconnected network centered on cellular response to DNA damage as modulated by CRREWs was highlighted (Figs. 2A, 3). This is consistent with a known effect of phthalate exposure resulting in DNA damage [50, 51].

Intriguingly, the fivefold paternally enriched CUL2, PHF10, EFCAB6, and FBRS, and twofold paternally enriched SMARCC1 and RARA are chromatin remodeler cofactors (Fig. 3, Additional file 5: Table S4) providing a foray into mechanism. As shown in Fig. 3, the majority of these eight DBP responsive and paternally provided RE-RNAs are involved in acetylation or deacetylation [27, 28, 38, 39, 44, 46, 52], except CUL2 and FBRS (Fig. 3). This highlights the importance of paternally derived acetylation factors during the final steps of spermatogenesis and, potentially, early embryogenesis. During spermatogenesis, acetylation of histone H4 is a critical step in replacing histones with protamine (reviewed in [53, 54]). Within 3 h of fertilization, the paternal chromatin will undergo transient hyperacetylation of histone H4 (reviewed in [53, 55]). To date, the molecular components integral to transient hyperacetylation remain elusive [53]. These RNAs may function in this transient hyperacetylation event.


Alterations in mouse sperm RNAs have been linked to offspring's metabolic health and stress response [2, 6, 8,9,10,11]. These studies have provided evidence in favor of the paternal origins of health and disease (POHaD) (reviewed in [12, 56, 57]). Recently, environmental exposures, including DBP and bisphenol A (BPA), and lifestyle factors such as BMI have been associated with alterations of epigenetic marks in sperm [4, 5, 12, 58,59,60,61], that is beginning to reconcile exposure and POHaD. Each of the CRREWs highlighted (CUL2, SMARCC1, PHF10, EFCAB6, FBRS, and RARA) alongside the non-CRREW CAMTA2 and PSME4 are paternally delivered as full-length RNAs ready for translation and early utilization in the fertilized oocyte. Perhaps these three CRREWs play a role directly following fertilization as the father’s chromatin is restructured or during syngamy. Interestingly, these genes are not represented within the human oocyte proteome [62], although CUL2, PSME4, and EFCAB6 are within the human sperm proteome [63]. As described above, they may encode early transient events like hyperacetylation in response to DBP exposure. On one hand, these RNAs are likely essential in functions prior to Embryonic Genome Activation, consistent with the MGI phenotype Ontology Annotations [64]. For example, CAMTA2, CUL2, PHF10, and SMARCC1 mouse knockdowns result in embryonic and/or preweaning lethality (Fig. 3, blue circles). On the other hand, PSME4 and RARA knockdowns impair male fertility due to several abnormalities related to spermatogenesis [64] (Fig. 3, blue circles). This emphasizes the importance of the sperm providing full-length transcripts and proteins at fertilization, as they may serve as the driving force behind the DBP-induced decreases in semen and embryo quality and subsequent increases in time to pregnancy.


DBP responsive REs

The differentially expressed REs summarized in Estill, MS et al. [5] from the crossover–crossback designed study were utilized (Additional file 2: Table S1). Men entered the study were on either a high-DBP-coated mesalamine at baseline (high-DBP study arm ( +), 112 semen samples, Fig. 1A) or non-DBP-coated mesalamine at baseline (background DBP study arm (−), 63 semen samples, Fig. 1A) [5]. It is important to note that both medications contain the same active pharmaceutical, mesalamine, and were exchangeably prescribed to IBD patients. They differed only in the presence of DBP in the coating. The 90-day intervals were designed to be reflective of a spermatogenic cycle and hence washout, when men would switch to the opposing drug from baseline to crossover, B1H (background DBP to high-DBP)/H1B (high-DBP to background DBP), then switched back from crossover to crossback, HB2 (high-DBP to background DBP)/BH2 (background DBP to high-DBP). Here, REs were evaluated as a function of this 90-day spermatogenic cycle and the duration of high-DBP exposure/withdrawal. This study was approved by the institutional review boards Partners Hospitals (Massachusetts General Hospital) protocol 2005P001631 and of Harvard T.H. Chan School of Public Health, Beth Israel Deaconess Medical Center, and Brigham and Women’s Hospital. The use of human tissue was approved by the Wayne State University Investigation Committee and carried out under the Wayne State University Human Investigation Committee IRB protocol 095701MP2E(5R).

Paternally provided REs

Paternally provided REs (generated from 7 non-IBD semen samples not exposed to high-DBP) were characterized from a total of 75,988 REs [11,386 RE–RE-containing RNAs (RE-RNAs)] identified within the zygote having a reads per kilobase per million (RPKM) > 10 [3]. As the zygote has yet to undergo embryonic genome activation, RE-RNAs present will be those provided directly by the sperm and oocyte [3]. Paternally enriched REs delivered at fertilization in which enrichment was at least fivefold higher when compared to the oocyte (fivefold paternally enriched) were described as having a median abundance > 25 RPKM in sperm, < 5 RPKM in the oocyte, and > 10 RPKM in the zygote [3]. This yielded a series of stringent REs that the father provides at fertilization.

To expand upon what may be paternally provided, an additional set of twofold paternally enriched REs were defined using a lower enrichment threshold for comparison (Fig. 1A, B). From the total 51,089 zygotic REs (10,277 RE-RNAs) independent of the paternally or maternally enriched REs previously defined [3], the paternal RE/maternal RE ratio was calculated in the following manner. If the RPKM of the zygotic RE was larger than the sum of the sperm and oocyte REs, the contribution of the sperm and oocyte equaled their respective RPKM abundance. If the zygote RPKM was less than the sum of the sperm plus oocyte REs, paternal contribution (Pc) was calculated as \(Pc = Z - \left( {\frac{Z}{{1 + \left( {\frac{s}{o}} \right)}}} \right)\), where Z represents the zygote REs RPKM, s represents the sperm RE RPKM and o represents the oocyte RE RPKM. Paternally contributed REs at a level twofold greater than the maternal contribution were termed twofold paternally enriched REs. While a series of these paternally provided REs are enriched in the sperm compared to the oocyte, some are specific to the sperm. For REs specific to the sperm, the RPKM in the oocyte was  < 2, the abundance value in which true RE presence cannot be confirmed. Each set of paternally provided REs was evaluated as a separate and combined RE list, as defined in Fig. 1.

Sample and transcript integrity

The transcript Integrity Index (TII) algorithm [26] was used to identify samples of similar quality [3, 5] using the 22 stable sperm-specific transcripts we previously defined [26]. The TII threshold was set at 50% of the transcript covered by at least 5 reads per million (RPM). Samples within the fourth quartile (Q4) were considered to have poor quality RNA [26]. Those samples passing TII were used to visually assess the paternally delivered RE corresponding RNAs of interest using the UCSC Genome Browser using Gencode version 36 [65]. RNAs were considered full-length if a minimum of 5 RPKM covered the transcript in all samples (7 paternally provided samples, 55 high-DBP study arm samples, 35 background DBP study arm samples). The 5 RPKM cutoff defines the minimum abundance in which there can be confidence in RE presence.

Gene ontology and statistical analysis

Enriched biological processes and pathways were evaluated using the Molecular Signature Database (MSigDB) version 7.5.1, employing the following collections: Hallmark gene sets (Hm), curated gene sets (C2), Gene Ontology (GO) gene sets (C5), and immunologic gene sets (C7). The collection C3: regulatory target gene set sub-category transcription factor (TF) targets were used to identify biologically corresponding TFs within the data. Each collection was considered separately. Thresholds were set to return the top 100 gene sets with a False Discovery Rate (FDR) q < 0.05 and a minimum of a two-gene overlap.

MSigDB analysis enables a maximum of 500 recognized genes per analysis. To evaluate the DBP responsive REs in Additional file 2: Table S1, each DBP comparison in Fig. 1A was separated into six groups based on the empirical (bootstrapped) p value that was generated using random resampling [5]. To group REs into the six empirical p value range groups, the REs within the largest DBP responsive comparison (all significantly associated REs within the(B1H comparison [5], Fig. 1B) were used. This would ensure that no comparison visualized in Fig. 1B would contain more than 500 unique RE-RNAs for MSigDB investigation. The six empirical p value range groups were as follows: group 1 = p < 0.013, group 2 = p between 0.013 and 0.023, group 3 = p between 0.023 and 0.032, group 4 = p between 0.032 and 0.041, group 5 = p between 0.041 and 0.045, and group 6 = p between 0.045 and 0.05. In the B1H, all correlated REs within Fig. 1B (3,651 [2,311 genes] REs), this segregated the REs as follows; group 1 = 711 REs (577 unique genes with 485 genes recognized), group 2 = 706 REs (566 unique genes with 486 genes recognized), group 3 = 718 REs (565 unique genes with 477 genes recognized), group 4 = 693 REs (551 unique genes with 470 genes recognized), group 5 = 451 REs (374 unique genes with 323 genes recognized), and group 6 = 372 REs (311 unique genes with 250 genes recognized).

EnrichR [66, 67], along with GeneCards ( [68] and the NIH Genetics Home Reference (, were utilized to assess gene function and disease associations. EnrichR categories of Pathways, Ontologies, and Diseases/Drugs were considered. Mediators of gene expression above TFs, considered Chromatin remodeler cofactors, RNA interactors, Readers, Erasers, and Writers (CRREWs). Briefly, CRREWs were identified from the curated list as described [4, 63], and key transcripts of interest were evaluated as part of the human sperm proteome.

The significance of proportional overlaps for TFs and CRREWs within the data was determined by the hypergeometric probability test with normal approximation from This provides a p value corresponding to a representation factor value indicating whether the overlap is significantly more or less than expected. A two-tailed t test for two samples of unequal variance was performed to calculate p values associated with the paternal/maternal contribution fold change using Microsoft 365 Excel (version 2202).

Availability of data and materials

Sequencing data for the DBP responsive REs are deposited as a GEO data set with the following accession number: GSE129216 and referenced in [5]. The REs utilized to determine the 5 × enriched and 2 × enriched paternally provided REs defined in [3]. Oocyte and zygote sequences were downloaded from the GEO database accessions GSE44183 and GSE71318.



Androgen receptor


Men transitioning from the high-DBP mesalamine to background DBP mesalamine (crossover visit)

B1 :

Men entering on the non-DBP mesalamine (background DBP baseline visit)


Background DBP (baseline visit) to high-DBP (crossover visit)

B2 :

Men returning to non-DBP mesalamine (background DBP crossback visit)

BH2 :

Background DBP (crossover visit) to high-DBP (crossback visit)

B1HB2 :

Background DBP study arm, baseline to crossover to crossback


Body mass index


Bisphenol A


Curated gene sets


Regulatory target gene set sub-category transcription factor targets


GO gene sets


Immunologic gene sets


Chromatin remodeler cofactors, RNA interactors, Readers, Erasers, and Writers




Environmental protection agency


False discovery rate


Gene ontology


Men transitioning from the non-DBP mesalamine to high-DBP mesalamine (crossover visit)

H1 :

Men entering on the high-DBP mesalamine (high-DBP baseline visit)


High-DBP (baseline visit) to background DBP (crossover visit)

H1BH2 :

High-DBP study arm, baseline to crossover to crossback

H2 :

Men returning to high-DBP mesalamine (high-DBP crossback visit)

HB2 :

High-DBP (crossover visit) to background DBP (crossback visit)


Histone deacetylase complexes


Hallmark gene sets


Irritable bowel disease


Monobutyl phthalate


Molecular signatures database

nBAF complex:

Neuron-specific chromatin remodeling complex


National Health and Nutrition Examination Survey

npBAF complex:

Neural progenitors-specific chromatin remodeling complex


RNA Element

Paternally provided REs:

5-Fold paternally enriched/twofold paternally enriched REs


Photodynamic therapy


Polycomb group


Paternal developmental origins of health and diseaseQ4: fourth quartile


RNA element discovery algorithm


RNA element


RE-containing RNA/gene


Reads per kilobase per million


Reads per million


Topological associated domain


Transcription factor


Transcript integrity index


  1. Ostermeier GC, Miller D, Huntriss JD, Diamond MP, Krawetz SA. Delivering spermatozoan RNA to the oocyte. Nature. 2004;429(6988):154.

    Article  CAS  Google Scholar 

  2. Godia M, Swanson G, Krawetz SA. A history of why fathers’ RNA matters. Biol Reprod. 2018;99(1):147–59.

    Article  Google Scholar 

  3. Estill MS, Hauser R, Krawetz SA. RNA element discovery from germ cell to blastocyst. Nucleic Acids Res. 2019;47(5):2263–75.

    Article  CAS  Google Scholar 

  4. Swanson G, Estill M, Diamond MP, Legro RS, Coutifaris C, Barnhart MD, et al. Human Chromatin remodeler cofactor interactor, eraser and writer sperm RNAs responding to obesity. Epigenetics. 2019.

    Article  Google Scholar 

  5. Estill M, Hauser R, Nassan FL, Moss A, Krawetz SA. The effects of di-butyl phthalate exposure from medications on human sperm RNA among men. Nat Sci Rep. 2019;9:1–3.

    CAS  Google Scholar 

  6. Gapp K, Jawaid A, Sarkies P, Bohacek J, Pelczar P, Prados J, et al. Implication of sperm RNAs in transgenerational inheritance of the effects of early trauma in mice. Nat Neurosci. 2014;17(5):667–9.

    Article  CAS  Google Scholar 

  7. Donkin I, Versteyhe S, Ingerslev LR, Qian K, Mechta M, Nordkap L, et al. Obesity and bariatric surgery drive epigenetic variation of spermatozoa in humans. Cell Metab. 2016;23(2):369–78.

    Article  CAS  Google Scholar 

  8. Grandjean V, Fourré S, De Abreu DAF, Derieppe M-A, Remy J-J, Rassoulzadegan M. RNA-mediated paternal heredity of diet-induced obesity and metabolic disorders. Sci Rep. 2015;5(1):1–9.

    Article  Google Scholar 

  9. de Castro BT, Ingerslev LR, Alm PS, Versteyhe S, Massart J, Rasmussen M, et al. High-fat diet reprograms the epigenome of rat spermatozoa and transgenerationally affects metabolism of the offspring. Mol Metab. 2016;5(3):184–97.

    Article  Google Scholar 

  10. Chen Q, Yan M, Cao Z, Li X, Zhang Y, Shi J, et al. Sperm tsRNAs contribute to intergenerational inheritance of an acquired metabolic disorder. Science. 2016;351(6271):397–400.

    Article  CAS  Google Scholar 

  11. Carone BR, Fauquier L, Habib N, Shea JM, Hart CE, Li R, et al. Paternally induced transgenerational environmental reprogramming of metabolic gene expression in mammals. Cell. 2010;143(7):1084–96.

    Article  CAS  Google Scholar 

  12. Marcho C, Oluwayiose OA, Pilsner JR. The preconception environment and sperm epigenetics. Andrology. 2020;8(4):924–42.

    Article  Google Scholar 

  13. Hauser R, Calafat AM. Phthalates and human health. Occup Environ Med. 2005;62(11):806–18.

    Article  CAS  Google Scholar 

  14. Lyche JL, Gutleb AC, Bergman Å, Eriksen GS, Murk AJ, Ropstad E, et al. Reproductive and developmental toxicity of phthalates. J Toxicol Environ Health, Part B. 2009;12(4):225–49.

    Article  CAS  Google Scholar 

  15. Nassan FL, Coull BA, Gaskins AJ, Williams MA, Skakkebaek NE, Ford JB, et al. Personal care product use in men and urinary concentrations of select phthalate metabolites and parabens: results from the environment and reproductive health (EARTH) study. Environ Health Perspect. 2017;125(8):087012.

    Article  Google Scholar 

  16. Hauser R, Duty S, Godfrey-Bailey L, Calafat AM. Medications as a source of human exposure to phthalates. Environ Health Perspect. 2004;112(6):751–3.

    Article  Google Scholar 

  17. Nassan FL, Korevaar TIM, Coull BA, Skakkebæk NE, Krawetz SA, Estill M, et al. Dibutyl-phthalate exposure from mesalamine medications and serum thyroid hormones in men. Int J Hyg Environ Health. 2019;222(1):101–10.

    Article  CAS  Google Scholar 

  18. Hernández-Díaz S, Mitchell AA, Kelley KE, Calafat AM, Hauser R. Medications as a potential source of exposure to phthalates in the US population. Environ Health Perspect. 2009;117(2):185–9.

    Article  Google Scholar 

  19. Bloom M, Whitcomb B, Chen Z, Ye A, Kannan K, Buck LG. Associations between urinary phthalate concentrations and semen quality parameters in a general population. Hum Reprod. 2015;30(11):2645–57.

    Article  CAS  Google Scholar 

  20. Wang Y-X, You L, Zeng Q, Sun Y, Huang Y-H, Wang C, et al. Phthalate exposure and human semen quality: Results from an infertility clinic in China. Environ Res. 2015;142:1–9.

    Article  CAS  Google Scholar 

  21. Wu H, Ashcraft L, Whitcomb BW, Rahil T, Tougias E, Sites CK, et al. Parental contributions to early embryo development: influences of urinary phthalate and phthalate alternatives among couples undergoing IVF treatment. Hum Reprod. 2017;32(1):65–75.

    CAS  Google Scholar 

  22. Louis GMB, Sundaram R, Schisterman EF, Sweeney A, Lynch CD, Kim S, et al. Semen quality and time to pregnancy: the longitudinal investigation of fertility and the environment study. Fertil Steril. 2014;101(2):453–62.

    Article  Google Scholar 

  23. Nassan FL, Coull BA, Skakkebaek NE, Williams MA, Dadd R, Mínguez-Alarcón L, et al. A crossover-crossback prospective study of dibutyl-phthalate exposure from mesalamine medications and semen quality in men with inflammatory bowel disease. Environ Int. 2016;95:120–30.

    Article  CAS  Google Scholar 

  24. Centers for Disease Control and Prevention, Fourth National Report on Human Exposure to Environmental Chemicals: CDC. 2021 Accessed 20 Aug 2022

  25. Shannon P, Markiel A, Ozier O, Baliga NS, Wang JT, Ramage D, et al. Cytoscape: a software environment for integrated models of biomolecular interaction networks. Genome Res. 2003;13(11):2498–504.

    Article  CAS  Google Scholar 

  26. Swanson GM, Estill MS, Krawetz SA. The transcript integrity index (TII) provides a standard measure of sperm RNA. Syst Biol Reprod Med. 2022;68(4):258–71.

    Article  CAS  Google Scholar 

  27. Niki T, Takahashi-Niki K, Taira T, Iguchi-Ariga SM, Ariga H. DJBP: a novel DJ-1-binding protein, negatively regulates the androgen receptor by recruiting histone deacetylase complex, and DJ-1 antagonizes this inhibition by abrogation of this complex. Mol Cancer Res. 2003;1(4):247–61.

    CAS  Google Scholar 

  28. Jin W. Novel insights into PARK7 (DJ-1), a potential anti-cancer therapeutic target, and implications for cancer progression. J Clin Med. 2020;9(5):1256.

    Article  CAS  Google Scholar 

  29. Medvedeva YA, Lennartsson A, Ehsani R, Kulakovskiy IV, Vorontsov IE, Panahandeh P, et al. EpiFactors: a comprehensive database of human epigenetic factors and complexes. Database (Oxford). 2015;2015:bav067.

    Article  Google Scholar 

  30. Fujiki R, Chikanishi T, Hashiba W, Ito H, Takada I, Roeder RG, et al. GlcNAcylation of a histone methyltransferase in retinoic-acid-induced granulopoiesis. Nature. 2009;459(7245):455–9.

    Article  CAS  Google Scholar 

  31. Schwartz YB, Pirrotta V. A new world of Polycombs: unexpected partnerships and emerging functions. Nat Rev Genet. 2013;14(12):853–64.

    Article  CAS  Google Scholar 

  32. Gao Z, Zhang J, Bonasio R, Strino F, Sawai A, Parisi F, et al. PCGF homologs, CBX proteins, and RYBP define functionally distinct PRC1 family complexes. Mol Cell. 2012;45(3):344–56.

    Article  CAS  Google Scholar 

  33. Banga SS, Peng L, Dasgupta T, Palejwala V, Ozer HL. PHF10 is required for cell proliferation in normal and SV40-immortalized human fibroblast cells. Cytogenet Genome Res. 2009;126(3):227–42.

    Article  CAS  Google Scholar 

  34. Staahl BT, Crabtree GR. Creating a neural specific chromatin landscape by npBAF and nBAF complexes. Curr Opin Neurobiol. 2013;23(6):903–13.

    Article  CAS  Google Scholar 

  35. Phelan ML, Sif S, Narlikar GJ, Kingston RE. Reconstitution of a core chromatin remodeling complex from SWI/SNF subunits. Mol Cell. 1999;3(2):247–53.

    Article  CAS  Google Scholar 

  36. Foster PM, Mylchreest E, Gaido KW, Sar M. Effects of phthalate esters on the developing reproductive tract of male rats. Hum Reprod Update. 2001;7(3):231–5.

    Article  CAS  Google Scholar 

  37. Izhar L, Adamson B, Ciccia A, Lewis J, Pontano-Vaites L, Leng Y, et al. A systematic analysis of factors localized to damaged chromatin reveals PARP-dependent recruitment of transcription factors. Cell Rep. 2015;11(9):1486–500.

    Article  CAS  Google Scholar 

  38. Song K, Backs J, McAnally J, Qi X, Gerard RD, Richardson JA, et al. The transcriptional coactivator CAMTA2 stimulates cardiac growth by opposing class II histone deacetylases. Cell. 2006;125(3):453–66.

    Article  CAS  Google Scholar 

  39. Martin M, Kettmann R, Dequiedt F. Class IIa histone deacetylases: regulating the regulators. Oncogene. 2007;26(37):5450–67.

    Article  CAS  Google Scholar 

  40. Qian M-X, Pang Y, Liu Cui H, Haratake K, Du B-Y, Ji D-Y, et al. Acetylation-mediated proteasomal degradation of core histones during DNA repair and spermatogenesis. Cell. 2013;153(5):1012–24.

    Article  CAS  Google Scholar 

  41. Ustrell V, Hoffman L, Pratt G, Rechsteiner M. PA200, a nuclear proteasome activator involved in DNA repair. EMBO J. 2002;21(13):3516–25.

    Article  CAS  Google Scholar 

  42. Srinivas H, Xia D, Moore NL, Uray IP, Kim H, Ma L, et al. Akt phosphorylates and suppresses the transactivation of retinoic acid receptor alpha. Biochem J. 2006;395(3):653–62.

    Article  CAS  Google Scholar 

  43. Li XS, Trojer P, Matsumura T, Treisman JE, Tanese N. Mammalian SWI/SNF-A subunit BAF250/ARID1 is an E3 ubiquitin ligase that targets histone H2B. Mol Cell Biol. 2010;30(7):1673–88.

    Article  CAS  Google Scholar 

  44. Hoyer J, Ekici Arif B, Endele S, Popp B, Zweier C, Wiesener A, et al. Haploinsufficiency of ARID1B, a member of the SWI/SNF-A chromatin-remodeling complex, is a frequent cause of intellectual disability. Am J Hum Genet. 2012;90(3):565–72.

    Article  CAS  Google Scholar 

  45. Pause A, Lee S, Worrell RA, Chen DY, Burgess WH, Linehan WM, et al. The von Hippel-Lindau tumor-suppressor gene product forms a stable complex with human CUL-2, a member of the Cdc53 family of proteins. Proc Natl Acad Sci USA. 1997;94(6):2156–61.

    Article  CAS  Google Scholar 

  46. Kamura T, Maenaka K, Kotoshiba S, Matsumoto M, Kohda D, Conaway RC, et al. VHL-box and SOCS-box domains determine binding specificity for Cul2-Rbx1 and Cul5-Rbx2 modules of ubiquitin ligases. Genes Dev. 2004;18(24):3055–65.

    Article  CAS  Google Scholar 

  47. Forrester NA, Sedgwick GG, Thomas A, Blackford AN, Speiseder T, Dobner T, et al. Serotype-specific inactivation of the cellular DNA damage response during adenovirus infection. J Virol. 2011;85(5):2201–11.

    Article  CAS  Google Scholar 

  48. Cukras S, Morffy N, Ohn T, Kee Y. Inactivating UBE2M impacts the DNA damage response and genome integrity involving multiple cullin ligases. PLoS ONE. 2014;9(7):e101844.

    Article  Google Scholar 

  49. Zhao W, Huang Y, Zhang J, Liu M, Ji H, Wang C, et al. Polycomb group RING finger proteins 3/5 activate transcription via an interaction with the pluripotency factor Tex10 in embryonic stem cells. J Biol Chem. 2017;292(52):21527–37.

    Article  CAS  Google Scholar 

  50. Duty SM, Singh NP, Silva MJ, Barr DB, Brock JW, Ryan L, et al. The relationship between environmental exposures to phthalates and DNA damage in human sperm using the neutral comet assay. Environ Health Perspect. 2003;111(9):1164–9.

    Article  CAS  Google Scholar 

  51. Hauser R, Meeker JD, Singh NP, Silva MJ, Ryan L, Duty S, et al. DNA damage in human sperm is related to urinary levels of phthalate monoester and oxidative metabolites. Hum Reprod. 2006;22(3):688–95.

    Article  Google Scholar 

  52. Heebøll S, Borre M, Ottosen PD, Andersen CL, Mansilla F, Dyrskjøt L, et al. SMARCC1 expression is upregulated in prostate cancer and positively correlated with tumour recurrence and dedifferentiation. Histol Histopathol. 2008;23(9):1069–76.

    Google Scholar 

  53. McLay DW, Clarke HJ. Remodelling the paternal chromatin at fertilization in mammals. Reproduction. 2003;125(5):625–33.

    Article  CAS  Google Scholar 

  54. Li J, Tsuprykov O, Yang X, Hocher B. Paternal programming of offspring cardiometabolic diseases in later life. J Hypertens. 2016;34(11):2111.

    Article  CAS  Google Scholar 

  55. Maalouf WE, Alberio R, Campbell KH. Differential acetylation of histone H4 lysine during development of in vitro fertilized, cloned and parthenogenetically activated bovine embryos. Epigenetics. 2008;3(4):199–209.

    Article  Google Scholar 

  56. Lacagnina S. The developmental origins of health and disease (DOHaD). Am J Lifestyle Med. 2020;14(1):47–50.

    Article  Google Scholar 

  57. Soubry A. POHaD: why we should study future fathers. Environ Epigenet. 2018;4(2):dvy007.

    Article  Google Scholar 

  58. Wu H, Estill MS, Shershebnev A, Suvorov A, Krawetz SA, Whitcomb BW, et al. Preconception urinary phthalate concentrations and sperm DNA methylation profiles among men undergoing IVF treatment: a cross-sectional study. Hum Reprod. 2017;32(11):2159–69.

    Article  CAS  Google Scholar 

  59. Tian M, Liu L, Zhang J, Huang Q, Shen H. Positive association of low-level environmental phthalate exposure with sperm motility was mediated by DNA methylation: a pilot study. Chemosphere. 2019;220:459–67.

    Article  CAS  Google Scholar 

  60. Zheng H, Zhou X, Li D-k, Yang F, Pan H, Li T, et al. Genome-wide alteration in DNA hydroxymethylation in the sperm from bisphenol A-exposed men. PLoS ONE. 2017;12(6):e0178535.

    Article  Google Scholar 

  61. Oluwayiose OA, Marcho C, Wu H, Houle E, Krawetz SA, Suvorov A, et al. Paternal preconception phthalate exposure alters sperm methylome and embryonic programming. Environ Int. 2021;155:106693.

    Article  CAS  Google Scholar 

  62. Virant-Klun I, Leicht S, Hughes C, Krijgsveld J. Identification of maturation-specific proteins by single-cell proteomics of human oocytes. Mol Cell Proteomics. 2016;15(8):2616–27.

    Article  CAS  Google Scholar 

  63. Castillo J, Jodar M, Oliva R. The contribution of human sperm proteins to the development and epigenome of the preimplantation embryo. Hum Reprod Update. 2018;24(5):535–55.

    Article  CAS  Google Scholar 

  64. Smith CL, Eppig JT. The mammalian phenotype ontology: enabling robust annotation and comparative analysis. Wiley Interdiscip Rev Syst Biol Med. 2009;1(3):390–9.

    Article  CAS  Google Scholar 

  65. Kent WJ, Sugnet CW, Furey TS, Roskin KM, Pringle TH, Zahler AM, et al. The human genome browser at UCSC. Genome Res. 2002;12(6):996–1006.

    Article  CAS  Google Scholar 

  66. Kuleshov MV, Jones MR, Rouillard AD, Fernandez NF, Duan Q, Wang Z, et al. Enrichr: a comprehensive gene set enrichment analysis web server 2016 update. Nucleic Acids Res. 2016;44(W1):W90–7.

    Article  CAS  Google Scholar 

  67. Chen EY, Tan CM, Kou Y, Duan Q, Wang Z, Meirelles GV, et al. Enrichr: interactive and collaborative HTML5 gene list enrichment analysis tool. BMC Bioinformatics. 2013;14:128.

    Article  Google Scholar 

  68. Stelzer G, Rosen N, Plaschkes I, Zimmerman S, Twik M, Fishilevich S, et al. The GeneCards suite: from gene data mining to disease genome sequence analyses. Curr Protoc Bioinformatics. 2016.

    Article  Google Scholar 

Download references


Not applicable.


Funding was from NIH Grants ES017285, ES009718 and ES000002 to RH and SAK. This work was further supported through the Charlotte B. Failing Professorship to SAK and the Robert J. Sokol, MD Chairship to JRP.

Author information

Authors and Affiliations



Planning and patient recruitment for the MARS (Mesalamine and Reproductive Health Study) study, which provided the DBP responsive REs, was by RH, FLN and JBF. RH and SAK designed and directed the original primary study. SAK directed the RNA isolation, characterization, and sequencing and data analysis. GMS performed data analysis of the presently described work. All authors read and approved the final manuscript.

Corresponding author

Correspondence to S. A. Krawetz.

Ethics declarations

Ethics approval and consent to participate

The semen samples were collected, processed, and analyzed at Wayne State University under the IRB protocols H-06-67-96 and HIC 095701MP2F.

Consent for publication

Not applicable.

Competing interests

The authors declare no competing financial interests but the following competing interest. FLN is a current employee and a shareholder of Biogen, Cambridge, MA. The original work on MARS Study (2 arm study) at Harvard T. H. Chan School of Public Health (HSPH), however, pre-dated the current employment. This manuscript does not mention any Biogen products or any of the disease states that Biogen is actively doing research in (to the coauthor’s knowledge). The authors declare no non-financial competing interests.

Additional information

Publisher's Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Supplementary Information

Additional file 1: Figure S1.

Chromatin remodeler cofactor, RNA interactor, reader, eraser and writer (CRREW) RNA Element (RE)-containing RNA (RE-RNA) visual integrity. Representative samples chosen for A) CUL2, B) SMARCC1 and C) PHF10. Integrity of DBP responsive and paternally provided CRREWs was determined using the UCSC Genome Browser Gencode version 41 track. Threshold for an RE-RNA to be considered full-length was set at a minimum of 5 Reads per Kilobase per Million (RPKM) across all transcript exons in all 7 paternally provided samples and all DBP responsive samples (high-DBP study arm (H1BH2), 55 samples; background-DBP study (B1HB2) arm: 35 samples).

Additional file 2: Table S1.

Number of RNA elements (REs) responsive to dibutyl-phthalate (DBP). REs were obtained from the publication Estill, MS et al. (2019b) (5). H1B; high-DBP (baseline visit) to background DBP (crossover visit), BH2; background DBP (crossover visit) to high-DBP (crossback visit), B1H; background DBP (baseline visit) to high-DBP (crossover visit), HB2; high-DBP (crossover visit) to background DBP (crossback visit).

Additional file 3: Table S2.

Enriched biological processes and pathways related to cellular stress, cell cycle, apoptosis, DNA damage response and gene regulation. RNA Element (RE) indicated in Fig. 1 panel B processes were used to query the Molecular Signature Database. Bolded text indicates the gene set is enriched in the paternally provided and DBP responsive RE-containing RNAs (RE-RNAs) evaluated. Italicized text indicates a gene set enriched upon DBP exposure addition and withdraw irrespective of original study arm. H1B; high-DBP (baseline visit) to background DBP (crossover visit), BH2; background DBP (crossover visit) to high-DBP (crossback visit), B1H; background DBP (baseline visit) to high-DBP (crossover visit), HB2; high-DBP (crossover visit) to background DBP (crossback visit), purple text; gene sets related to cellular stress, brown text; gene sets related to the cell cycle, green text; gene sets related to apoptosis and cell death, red text; gene sets related to DNA damage response (DNA damage response to UV exposure), blue text; gene sets related to gene regulation., indicates no enrichment in the gene set.

Additional file 4: Table S3.

Number of dibutyl-phthalate (DBP) responsive transcription factor (TF) binding site gene sets and Chromatin remodeler cofactors, RNA interactors, Readers, Erasers and Writers (CRREWs). A) Enriched TF binding site gene sets. Enriched gene sets were separated based on the identification of having a known TF reported to bind. B) Number of DBP responsive CRREWs. Gene lists were generated from those RNA Element (RE)-containing RNAs (RE-RNAs) represented within Fig. 1. Representation factor and p value were determined by hypergeometric probability test.

Additional file 5: Table S4.

All Di-butyl phthalate (DBP) responsive and paternally provided CRREWs within sperm. H1B; high-DBP (baseline visit) to background DBP (crossover visit), BH2; background DBP (crossover visit) to high-DBP (crossback visit), B1H; background DBP (baseline visit) to high-DBP (crossover visit), HB2; high-DBP (crossover visit) to background DBP (crossback visit), paternally provided, RE-containing RNAs (RE-RNAs) that are fivefold paternally enriched or twofold paternally enriched.

Additional file 6: Table S5.

REs unique to the zygote that are twofold paternally enriched. A) DBP responsive and 2x paternally enriched. B) 2x paternally enriched but not DBP responsive. REs required the paternal contribution to be > twofold the maternal, or the maternal RE abundance to be < 2 RPKM and paternal RE abundance < 25 and > 2 RPKM. #DIV/0 indicates contribution is solely from the father. Green fill indicates the RE is shared between two genes in the sperm contributed set.

Additional file 7: Table S6.

Number of paternally provided Chromatin remodeler cofactors, RNA interactors, readers, erasers and writers (CRREWs). A) Number of paternally provided CRREWs. B) Number of dibutyl-phthalate (DBP) responsive, paternally provided CRREWs. Representation factor and p value were determined by hypergeometric probability test. Paternally provided CRREWs include those that are fivefold paternally enriched and twofold paternally enriched.

Additional file 8: Table S7.

Integrity of specifically paternally provided and di-butyl phthalate (DBP) responsive transcripts. RNA Element (RE)-containing RNAs (RE-RNAs) represented respond to DBP addition and subtraction in the same direction. Visual inspection is based on the UCSC Genome Browser. Bold text indicates the RE is shared between the paternally provided and DBP responsive samples, blue text indicates the DBP responsive RE was associated with the baseline to crossover (B1H or H1B) comparison and red text indicates the DBP responsive RE was associated with the crossover to crossback comparison (HB2 or BH2). Paternally provided indicates samples in which the fivefold paternally enriched and twofold paternally enriched REs were obtained. B1HB2; background DBP study arm baseline to crossover to crossback, H1BH2; high-DBP study arm baseline to crossover to crossback. Mean RE abundance is in Reads per Kilobase per Million (RPKM).

Additional file 9: Table S8.

Gene Ontology of dibutyl-phthalate (DBP) responsive, paternally provided RNA Element (RE)-containing RNAs (RE-RNAs) not within the enriched biological processes. Ontology was assigned using EnrichR. Bold text; indicates function related to the enriched biological processes highlighted in Fig. 2B, italics; indicate an ontology that is shared by at least two RE-RNAs, na; no ontology reported, not within EnrichR; the RE-RNA is not recognized by EnrichR, …; no relevance to enriched biological processes highlighted in Fig. 2B, B1HB2; background DBP study arm baseline to crossover to crossback, H1BH2; high-DBP study arm baseline to crossover to crossback, Paternally provided; fivefold paternally enriched and twofold paternally enriched RE-RNAs.

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit The Creative Commons Public Domain Dedication waiver ( applies to the data made available in this article, unless otherwise stated in a credit line to the data.

Reprints and Permissions

About this article

Verify currency and authenticity via CrossMark

Cite this article

Swanson, G.M., Nassan, F.L., Ford, J.B. et al. Phthalates impact on the epigenetic factors contributed specifically by the father at fertilization. Epigenetics & Chromatin 16, 3 (2023).

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI:


  • Phthalates
  • Sperm RNA
  • Paternal contribution
  • Chromatin modifiers