Genetic resistance to DEHP-induced transgenerational endocrine disruption

Di(2-ethylhexyl)phthalate (DEHP) interferes with sex hormones signaling pathways (SHP). C57BL/6J mice prenatally exposed to 300 mg/kg/day DEHP develop a testicular dysgenesis syndrome (TDS) at adulthood, but similarly-exposed FVB/N mice are not affected. Here we aim to understand the reasons behind this drastic difference that should depend on the genome of the strain. In both backgrounds, pregnant female mice received per os either DEHP or corn oil vehicle and the male filiations were examined. Computer-assisted sperm analysis showed a DEHP-induced decreased sperm count and velocities in C57BL/6J. Sperm RNA sequencing experiments resulted in the identification of the 62 most differentially expressed RNAs. These RNAs, mainly regulated by hormones, produced strain-specific transcriptional responses to prenatal exposure to DEHP; a pool of RNAs was increased in FVB, another pool of RNAs was decreased in C57BL/6J. In FVB/N, analysis of non-synonymous single nucleotide polymorphisms (SNP) impacting SHP identified rs387782768 and rs29315913 respectively associated with absence of the Forkhead Box A3 (Foxa3) RNA and increased expression of estrogen receptor 1 variant 4 (NM_001302533) RNA. Analysis of the role of SNPs modifying SHP binding sites in function of strain-specific responses to DEHP revealed a DEHP-resistance allele in FVB/N containing an additional FOXA1-3 binding site at rs30973633 and four DEHP-induced beta-defensins (Defb42, Defb30, Defb47 and Defb48). A DEHP-susceptibility allele in C57BL/6J contained five SNPs (rs28279710, rs32977910, rs46648903, rs46677594 and rs48287999) affecting SHP and six genes (Svs2, Svs3b, Svs4, Svs3a, Svs6 and Svs5) epigenetically silenced by DEHP. Finally, targeted experiments confirmed increased methylation in the Svs3ab promoter with decreased SEMG2 persisting across generations, providing a molecular explanation for the transgenerational sperm velocity decrease found in C57BL/6J after DEHP exposure. We conclude that the existence of SNP-dependent mechanisms in FVB/N inbred mice may confer resistance to transgenerational endocrine disruption.

Introduction Di-(2-ethylhexyl) phthalate (DEHP; CAS No. 117-81-7) is a reproductive toxicant and an endocrine disruptor (ED) ubiquitously found in the environment. Accumulated data demonstrate that DEHP interferes with sex steroid hormone signaling pathways (SHP). DEHP and its principal metabolite named mono-(2-ethylhexyl) phthalate (MEHP; CAS No. 4376-20-9) decrease the testosterone produced by testes and interact at the molecular level with the androgen (AR), estrogen (ER) and peroxisome proliferator-activated receptors (PPARs) [1,2]. Prenatal exposure to DEHP causes androgen deficiency during embryogenesis in both animals and humans [3,4]. The anogenital distance (AGD), a marker of fetal androgen exposure [5], was shortened in boys born from DEHP-exposed mothers and was reduced in rodents prenatally exposed to DEHP [6][7][8]. Therefore, the long-term toxicological impacts of prenatal exposure to DEHP are of high concern.
We injected per os 300 mg/kg/day DEHP to pregnant mice during embryonic (E) days (E9-19), and measured male fertility parameters at adulthood. The dose was chosen from a previous study and appears to be relevant for extreme human exposure. In fact, the dose of DEHP effectively reaching the mice fetus in the present study was estimated at 190 μg/kg/day and is comparable with the 233 μg/kg/day of median daily intake of DEHP in neonates treated in intensive care units [9]. First, 55% of ingested DEHP is absorbed, whereas DEHP and its derivatives are predominately excreted in the urine. In addition, approximately 20-25% of absorbed DEHP cannot pass the gastrointestinal tract barrier of the pregnant animal or mother, and is excreted in the feces (ToxGuide for DEHP). Thus, a fraction of excreted DEHP is not able to reach the embryos in pregnant females. In fact, only 0.03% of the initial dose of 14 C-labelled DEHP, orally administrated to pregnant mice at 8 days of gestation, was recovered in the fetuses when monitoring radioactivity [10]. Among the 9 mg of DEHP that were given per pregnant mice per days, the reconstructed dose of DEHP effectively received by the fetus is estimated at 190 μg/kg/day; 0.27 μg of the initial dose reaches the fetal tissue weighting 1.4 � 10-3 kg. That dose is lower than the median daily intake of DEHP calculated in infants in the high-intensiveness product use group. This dose was estimated to range from 233 to 352 μg/ kg/day based on MEHHP and MEOHP concentrations recovered in the urines of the preterm infants exposed to DEHP-containing medical products [9]. However, the metabolites that reach the embryos may differ, with DEHP metabolites produced by the exposed mother on one hand, and direct leaching of DEHP from the medical products in the blood circulation of the neonates on the other hand. As a result, a decreased sperm count was observed in the C57BL/6J strain, but not in FVB/N mice, indicating that the latter seem to be resistant and the former sensitive to DEHP [11]. Previously, heterogeneity explained by strains was reported in DEHP-exposed mice [6]. We think that resistance to prenatal exposure to DEHP may imply genetic variations affecting the direct or indirect targets of DEHP, in enzymes responsible for excretion of DEHP, or in DNA sequences recognized by the hormones that are affected by DEHP. The exposure mechanism implies that DEHP orally injected in the mouse mouth cavity passes into the digestive tract of the pregnant female, is metabolized notably in MEHP, enters the blood circulation of pregnant mice, reach the embryos through the cord blood, and exert anti-androgenic and pro-estrogenic activities. The main suspected impact in the male fetus is a decreased production of testosterone. This may trigger long-term impacts on the future sexual health at adulthood, such as in the testicular dysgenesis syndrome (TDS) [12,13].
The exposure time frame extends from the primordial germ cell (PGC) migration period (~E10.5) and covers the gonadal differentiation period initiated between E11 and E12 in the developing mice embryos, representing a susceptible time window for androgen interference [14]. PGCs are characterized by a specific pattern of DNA methylation and transcription controlling the cell fate in the lineage [14]. PGCs differentiate themselves in male germ cells during the fetal androgen-dependent sexual differentiation of the somatic gonads to testis. At puberty, a peak of testosterone triggers the spermatogenesis process starting from the spermatogonia. The spermatogenesis involves transcriptional patterns specific to the successive developmental stages [15]. Thus, endocrine disruption taking place during PGCs migration may probably alter DNA methylation and RNA expression status in mature sperm at adulthood.
It has previously been reported that aberrant DNA methylation or dysregulated RNAs in the sperm of adult mice were able to support transgenerational inheritance [16,17]. Therefore, we studied the transmission of the biological alterations to the subsequent generations taking into account the epigenetic status of genes encoding the dysregulated sperm RNAs. Active transcription is believed to be absent in mature spermatozoa, but RNAs remain present deeply embedded in the sperm nucleus [18]. At least part of these paternal RNAs is transmitted to the egg at fertilization, which may help to understand "unexplained male-factor infertility" [19]. Moreover, a functional AR is expressed in both X and Y carrier spermatozoa and dihydrotestosterone (DHT) is present in the seminal fluid [20]. Additionally, steroid receptors and their ligands may impact on male gamete functions [21].
We postulate that FVB/N resistance or alternatively C57BL/6J susceptibility, may involve genomic variations affecting the androgen signaling occurring in a strain-specific manner. Both DNA methylation and RNA expression may be affected in the sperm of mice prenatally exposed to DEHP, whereas SNP may influence both DNA methylation and gene expression [22]. Comparison between both genomes performed by Wong K. et al were incorporated in the present study [23]. Genomes comparison revealed that the AR is not directly affected by polymorphisms in FVB/N compared with the C57BL/6J reference genome. However, the nonsynonymous single nucleotide polymorphism (SNP) rs29315913 ("A" in FVB/N and "C" in C57BL/6J) was found in one variant of ESR1 and rs387782768 ("C" in FVB/N and "T" in C57BL/6J) was found in the Forkhead box A3 (FOXA3) transcription factor required for testicular steroidogenesis [23,24]. Both factors are important for male fertility and are key player of SHP. Briefly, ESR1 is involved in the regulation of male reproductive organs, notably the efferent ductus and the epididymis, and its loss results in impaired ion transport and water reabsorption, with production of abnormal sperm characterized by abnormal flagellar coiling and increased incidence of spontaneous acrosome reactions [25,26]. The Forkhead box A genes (FOXA1, FOXA2 and FOXA3) encode pioneer transcription factors, i.e. the first detectable factors to engage target sites in chromatin, thus facilitating the further binding of the AR in the nucleus to the targeted recognized DNA sequences [27].
The principal aim of this genome-environment interaction study focalized on male germ cells was to identify SNP varying between FVB/N and C57BL/6J genomes and responsible for the different phenotypes. We postulate that some SNP interacting with DEHP exposure may change the RNA content and the DNA methylation status in the sperm in a strain-specific manner with putative transgenerational inheritance.

Phenotypic impact of in utero exposure to DEHP in C57BL/6J and FVB/N mice
The phenotypic changes affecting male fertility parameters induced by prenatal exposure to DEHP confirm the FVB/N resistance found in our earlier study using an independent method [11] and using data partially previously obtained in C57BL/6J [13]. In FVB/N, AGD, testes weight, sperm concentration and curvilinear velocity (VSL) were not affected by prenatal exposure to DEHP and remain stable also in the second generations of mice generated from the prenatally exposed F1 males (Fig 1A and 1B). The only exception concerns sperm velocities with decreased average path velocity (VAP) and VSL in FVB.D300.F1 compared to FVB.CTL. F1. A return to control values was observed in FVB.D300.F2 (Fig 1B). The experiments were stopped at F2 in FVB/N due to the absence of any remaining impact on all tested parameters.
Overall, the results demonstrated a strain-specific difference in the impact of prenatal exposure to DEHP in FVB/N compared with C57BL/6J inbred mice. An inbred mice strain is a population of mice clones at the genetic level resulting from a process of at least 20 sequential generations of brother-sister mating. They also strongly suggest the persistence of a detrimental impact across generations in the case of DEHP-susceptibility and revealed a surprising deterioration of sperm velocity apparently inherited in a transgenerational manner (from C57. F1.D300 to C57.F3.D300).

Sperm transcriptome variations induced by prenatal exposure to DEHP in C57BL/6J and FVB/N mice
SHP were shown to be affected by DEHP and its metabolites in previous reports. Both androgens and estrogens exerted their biological effects on gene transcription upon binding to their respective receptors (AR and ESR1). Here, we postulate that prenatal exposure to DEHP may dysregulate sperm RNAs in a persistent manner. Total RNA was extracted from sperm samples analyzed by CASA and the RNA content was analyzed using all-RNA-seq (Fig 2A and S1  Table). One differential analysis was performed using multiple pairwise comparisons between the conditions C57.CTL.F1, C57.D300.F1, FVB.CTL.F1 and FVB.D300.F1 in order to take into account the impact of the strains on sperm RNA content in the absence or presence of DEHP exposure. Additionally, as almost all parameters tested remained affected at F2 in C57BL/6J (Fig 1), we also included the C57.D300.F2 condition.
This approach resulted in genome � environment � generation interactions transcriptomic analysis identifying 62 RNAs as the most differentially expressed RNA among the different conditions ( Fig 2B). Importantly, three different patterns of expression changes were found (Fig 2C, 2D and 2E). The first pattern showed increased RNA levels specifically in FVB.D300. F1 compared with FVB.CTL.F1. This pattern presenting an FVB/N-specific and DEHP-mediated induction of RNA expression was associated with FVB/N resistance (Fig 2C). The second pattern was associated with generational changes occurring between the first and second generations in the C57BL/6J exposed lineage (Fig 2D). The third pattern revealed decreased levels of various RNA specifically in C57.D300.F1 compared with C57.CTL.F1, without changes in FVB/N. The latter pattern was associated with C57BL/6J susceptibility to DEHP (Fig 2E). Overall, the 62 differentially expressed RNA across the tested conditions were mainly regulated by hormones in sperm or in male reproductive tissues ( Fig 2F and Table 1). The vast majority (80%) present increased expression during spermatogenesis ( Fig 2G); androgen was the most represented regulatory hormone ( Fig 2H). In addition, pathways highly relevant to sperm physiological functions and DEHP toxicity were affected across the tested conditions (S2 Table).
Strain-specific non-synonymous SNPs affecting AR, ESR1 and FOXA1-3 sperm RNA levels SNP variations between FVB/N and C57BL/6J affecting SHP may explain the strain-specificity of the DEHP-induced alterations in sperm RNAs (Fig 2). Non-synonymous SNP rs387782768 in FOXA3 and rs29315913 in ESR1 variant 4 genes were characterized previously [23]. No other relevant SNP affecting SHP could be directly found in both strains. Transcriptomic results of this study revealed that both SNPs were associated with strainspecific differences of particular RNA expression levels, independent of DEHP exposure ( Fig 3A). First, rs387782768 was associated with the absence of FOXA3 sperm transcripts in FVB/N. Second, rs29315913 was associated with increased levels of ESR1 variant 4 sperm transcript in FVB/N. Other tested transcription factors not affected by non-synonymous SNP (FOXA1, FOXA2, AR and other ESR1 variants) did not differ in their expression levels between both strains. An upregulation of AR expression was also observed in C57.D300.F2 that did not depend directly on a SNP, apparently reflecting a generational impact in the exposed lineage of the C57BL/6J strain.

Effects of DEHP on promoter methylation and on gene expression in the sperm of both strains
Changes in RNA levels observed in sperm upon prenatal exposure to DEHP may involve alteration in epigenetic mechanisms. Among these, aberrant DNA methylation(s) in the male germ cell lineage caused by prenatal DEHP exposure may result in persistent alterations of sperm RNA patterns detectable at adulthood. To test this hypothesis, all-RNA-seq data were merged to previously acquired Methyl-CpG-binding domain sequencing (MBD-seq) data for all conditions tested at F1 also performed in quintuplicates (C57.CTL.F1, C57.D300.F1, FVB.CTL.F1 and FVB. D300.F1) [11]. The first observation is that transcript abundance was negatively correlated with promoter methylation in the sperm, consistent with a repressive impact of promoter-methylation on transcription named silencing ( Fig 3B) [112]. Promoter methylation was then tested specifically in the RNA that produced both patterns associated with either DEHP-susceptibility or DEHP resistance (Fig 2C and 2E and Table 1). The targets for DEHP-susceptibility, characterized by DEHP-induced decreased RNA expression levels, were all associated with a DEHPinduced increased promoter methylation in C57BL/6J (Fig 3B). These combined results are consistent with a DEHP-induced epigenetic silencing of the DEHP-susceptible targets identified only in C57BL/6J, as the same targets were not affected by DEHP in FVB/N, neither in their expression, nor in their promoter methylations levels ( Fig 3B). The targets for DEHP-resistance, upregulated by DEHP specifically in FVB/N, were also associated, paradoxically, with promoter methylation increases ( Fig 3B). The targets for DEHP-resistance were not affected by DEHP in C57BL/6J, neither in their expressions, nor in their promoter methylation levels ( Fig 3B).
We postulated that the decreased expression of androgen-regulated genes in particular may be directly explained by the anti-androgenic activity of DEHP in C57BL/6J, whereas the increased expression of genes in FVB/N may involve an additional pro-estrogenic activity of DEHP. Therefore, the binding sites for both estrogen and androgen signaling systems were analyzed in the promoters of the identified targets for both DEHP-susceptibility and DEHPresistance ( Fig 3C). Results revealed that receptor binding sites for androgen signaling were higher than for estrogen signaling in the promoters of targets for DEHP-susceptibility ( Fig  3C), whereas equal numbers of ER and AR binding sites were detected in the promoters of targets for DEHP-resistance ( Fig 3C). These results are compatible with anti-androgenic silencing and pro-estrogenic induction of the strain-specific targets.
Overall, these results demonstrate increased promoter methylations in strain-specific affected targets for DEHP-susceptibility in C57BL/6J, and for DEHP-resistance in FVB/N, with differences in the number of ER over AR binding sites found in the promoters of these strain-specific targets.

Impact of SNPs on strain-specific DEHP-induced sperm RNA changes
SNP variations affecting binding sites for FOXA1-3, AR and ESR1 and transcription factors were investigated in function of the DEHP-mediated, strain-specific sperm RNA changes ( Fig  4). Such an approach may allow the identification of alleles associated with either DEHP- Rprl3

DEHPresistance
All cells NA/NA/NA tRNA maturation - [110] Overall, these results demonstrate that prenatal exposure to DEHP may alter the male meiotic program in a strain-specific manner, producing distinct patterns of expression changes and mainly implying dysregulation in androgen-regulated sperm RNAs. https://doi.org/10.1371/journal.pone.0208371.t001 Transgenerational endocrine disruption susceptibility or DEHP-resistance if several SNPs are co-localized or if various co-localized genes are affected. A total of 6 SNPs showed the highest combined impacts on both the strain-specific transcriptional changes induced by prenatal exposure to DEHP measured with Z values, as well as on SHP-motifs specificities measured with δ values (Fig 4A).
We were able to identify rs30973633 as affecting a binding motif for all FOXA1-3 transcription factors in that this motif is present in the FVB/N and absent in the C57BL/6J allele. Moreover, the FVB/N-specific binding motif is associated via a probable DEHP-mediated induction of FOXA1-3 ( Fig 3A), with the activation of various beta-defensins localized around rs30973633 ( Fig 4B and Table 2). Precisely, prenatal exposure to DEHP induced increased expressions in Defb42, Defb30, Defb47 and Defb48, but not Defb43. None of these genes were affected in their expression after prenatal exposure to DEHP in C57BL/6J (Fig 4B). This DNA region displaying an FVB/N specific induction of RNA expression combined with an FVB/N specific binding site for FOXA1-3 was defined as an "FVB/N resistance allele" to prenatal DEHP exposure.
The second result was the identification of 5 SNP-dependent binding motifs co-localized in the same DNA region encoding the Svs genes ( Fig 4C). All these SNP-dependent binding motifs were specific to the C57BL/6J allele. More precisely, in the latter, an AR binding site in Svs4 was affected by two different SNPs (rs28279710 and rs46648903), another AR binding site between Svs3a and Svs6 was affected by rs46677594, an ESR1 binding site upstream of Svs3b was affected by rs32977910, and a FOXA2 binding site between Svs3a and Svs6 was affected by rs48287999. All these SNP-dependent binding sites were thus specific of the C57BL/6J allele and all the corresponding genes were silenced by DEHP in C57BL/6J only (Fig 4C and  Table 2). This region was defined as the "C57BL/6J susceptibility allele" to prenatal exposure to DEHP.
To summarize, the genome-wide analysis performed resulted in the characterization of the FVB/N resistance allele, involving 4 out of 5 beta-defensins, with all induced by DEHP in FVB/N only, as well as a C57BL/6J susceptibility allele encoding 6 Svs genes that were all silenced by DEHP in C57BL/6J only. Correlations were assessed using Kendall's tau and Spearman's rho rank correlation tests. Blue points represent the genes showing the pattern of DEHP-resistance involving Gpx5, Spag11b, Defb30, Lcn8, Lcn9, Lcn5, Spint4 and Rprl3, that was absent in the MBD-seq, but without 9230104L09Rik. Red points represent the genes showing the pattern of DEHP-susceptibility involving Sva, Svs2, Svs3a, Svs3b, Svs4, Svs5, Svs6, Spink1 and Saa2, without Pate4 and 950003J23Rik, both absent in MBD-seq. (C) Binding sites for AR (green), ER (pink) and FOXA1-3 transcription factors (black, dark gray and white) detected in promoters of genes associated with DEHPsusceptibility (red names) and DEHP-resistance (blue names). The number of AR (n = 6) binding sites were higher than ER (n = 2) sites in promoters of genes (n = 8) silenced by DEHP in C57BL/6J and associated with susceptibility. Equal numbers of ER (n = 6) and AR (n = 6) binding sites were detected in the genes (n = 9) induced by DEHP in FVB/N and associated with resistance. Binding sites are shown as triangles orientated differently depending on the DNA strands and with thickness reflecting the score. Sequences were extracted from -2000 to 0 bp with gene names from Mus musculus GRCm38. Analysis was performed using the matrix-scan program of Rsat [111]. https://doi.org/10.1371/journal.pone.0208371.g003 Transgenerational endocrine disruption

CpG methylation analysis in C57BL/6J-susceptibility allele combined with quantifications of the semenogelins (SEMG) production in sperm
Two independent targeted methods were applied to validate the previous findings on the susceptibility allele of C57BL/6J (Fig 5A). First, measures of CpG methylation levels were performed using bisulfite pyrosequencing in the promoters of both Svs2 and Svs3ab ( Fig 5B). Second, expression levels of SEMG1 and SEMG2 proteins encoded by Svs2 and Svs3ab, respectively, were analyzed by Western blots (Fig 5C and 5E). DNA and protein samples required for these analyzes were extracted by sequential precipitations from the TRIZOL interphase and organic phase that remain in frozen sperm samples after RNA extraction.
The bisulfite pyrosequencing results demonstrated increased methylation levels in three of four CpG sites localized in promoters of both Svs2 and Svs3ab in each of the exposed lineages affecting different generations (Svs2-CpG 1 methylation increased in C57.D300.F1 and in FVB. D300.F2; Svs3ab-CpG 1 methylation increased in both C57.D300.F3 and in FVB.D300.F2; Svs3ab-CpG 2 methylation increased in FVB.D300.F2), whereas Svs2 CpG 2 was always fully methylated ( Fig 5B). Methylation levels in both CpG sites of Svs3ab promoter were also increased gradually across the exposed lineage in C57BL/6J as revealed by correlations analyses (Fig 5B). Electrophoresis of total protein extracts combined with Western bolt revealed that ). δ = score C57BL/6J -score FVB/N . The corresponding logos above each graph represent the tested motif. Dashed lines represent thresholds. (B) FVB/N-specific DEHP-induced increased expression in the β-defensin loci associated with an FVB/N-specific FOXA1-3 binding motif due to rs30973633. Upon prenatal exposure to DEHP, increased expression levels are recorded for Defb42, Defb30, Defb47 and Defb48, but not for Defb43 in FVB/N, without changes in C57BL/6J. (C) C57BL/6J-specifc DEHP-induced expressional changes in the Svs loci associated with C57BL/6J-specific binding sites for AR, ESR1 and FOXA2 due to SNPs rs28279710, rs46648903, rs46677594, rs32977910 and rs46677594. In utero exposure to DEHP resulted in decreased expression across all Svs genes in C57BL/6J with no change in FVB/N. (BC) Delta values in parentheses represent score differences between both alleles for the respective motif, taking the SNP into account. Stars represent statistically significant differences of expression levels, according to the non-parametric pairwise Wilcoxon test with Benjamini & Hochberg correction for multiple testing. In part C, statistically significant differences between both strains are not shown.
https://doi.org/10.1371/journal.pone.0208371.g004 SEMG2 is the most expressed proteins in sperm samples, whereas SEMG1 was detected at very low levels ( Fig 5C). Western blot quantifications revealed decreased SEMG2 specifically in the exposed linage of C57BL/6J at F1 and F3, without changes in FVB/N ( Fig 5D). SEMG1 levels were unchanged across the different conditions tested in both backgrounds ( Fig 5E).
Overall, the DEHP-induced increased methylation levels recorded at the promoters of Svs3ab in the exposed lineage of the strain C57BL/6J were compatible with the observed decreases of SEMG2 expression in C57.D300.F1 and C57.D300.F3 conditions (Fig 5D). In addition, the results provide a molecular explanation for the DEHP-induced transgenerational decrease in sperm velocities recorded in the C57BL/6J lineages (Fig 1D, VAP, VCL and VSL). The latter might be due to the transgenerational inheritance of decreased SEMG2 production mediated by increased promoter methylation in Svs3ab. Indeed, it is recognized that a fundamental physiological function of semenogelins is the control of sperm motility [113]. DEHP effects on Svs3ab in C57BL/6J strain are coherent with the general consensus, showing an increased methylation, a decreased RNA transcription, and a decreased level of the corresponding protein, SEMG2. DEHP effects on Svs2 in this same strain are also coherent, showing an increased methylation in the promoter of Svs2 and a decreased transcription of Svs2 gene. The absence of effect on the corresponding protein, SEMG1, might be due to technical difficulties associated with the low level of this protein. The results obtained in the FVB/N strain were also coherent with the lack of effect of DEHP, except for the increase in both promoter methylation observed in F2 offspring.

Discussion
C57BL/6J mice prenatally exposed to DEHP present anti-androgenic and pro-estrogenic symptoms, whereas FVB/N mice are phenotypically clearly unaffected (Fig 1). Both SNP affecting directly SHP (rs29315913 and rs387782768) were found to be functional (Fig 3). First, rs29315913 affects the ligand-binding domain of one of the five ESR1 proteins by replacing glycine 447 for a valine due to a C to A conversion in FVB/N [23]. Consequently, NM_001302533 RNA is more expressed in FVB/N. Second, rs387782768 is associated with the absence of FOXA3 in FVB/N (Fig 3A).

FVB/N resistance to DEHP
In FVB/N, a pro-estrogenic impact of DEHP mediated by ESR1 variant 4 (NM_001302533) was detected. NM_001302533 was more expressed in FVB/N containing the SNP-affected Transgenerational endocrine disruption ligand binding domain ( Fig 3A) and a 3-fold higher number of ESR1 binding sites was detected across promoters of FVB/N-specific targets (Fig 3C). The involvement of ER in DEHP-mediated toxicity has been reported previously [114]. The pool of targets associated with resistance to DEHP in FVB/N were previously involved in sperm resistance mechanisms (Table 1), against bacterial infections or innate host defense (Defb30, Spag11b), in protective mechanisms against oxidative stress (Gpx5), and in the clearance of exogenous compounds (Lcn5, Lcn9, Lcn5) [44,115] in accordance with previously reports [116,117]. We identified a FVB/N-specific binding motif associated with a DEHP-induced induction in loci encoding several β-defensins (Fig 4B). Precisely, rs30973633 was shown to add a FOXA1-3 putative binding motif in the FVB/N allele encoding various beta-defensins overexpressed by DEHP specifically in that strain, with a positive impact of DEHP exposure on FOXA1-3 transcript levels. β-defensins have been shown to play a dual role in the sperm, protecting the sperm from bacterial infections, but also acting in a positive manner on sperm motility [118][119][120][121][122][123][124].

C57BL/6J susceptibility to DEHP
In C57BL/6J, the prenatal exposure to DEHP affected all CASA parameters (Fig 1D), reflecting perturbations of SHP such as for the decreased AGD [5], in consistency with meta-analyses [6]. Decreased sperm counts are most probably a consequence of alterations in the proliferation of the fetal male germ cell precursors in the gonadal differentiation period at the time of exposure, or decreased proliferation of Sertoli cells [125]. Decreased testis weight may be due to an alteration in the initial proliferation of all or any gonadal cells lineages or may also involve a long-term pro-estrogenic impact [126]. Interestingly, ER knock-out mice were infertile and showed increased fluid pressure within testes due to the accumulation of the fluid at its production site with decreased sperm counts [127]. Overall, DEHP-mediated anti-androgenic and pro-estrogenic symptoms that remain poorly understood were detected in C57BL/ 6J, but absent in FVB/N, consistent with an intrinsic resistance of the strain to this recognized endocrine disruptor (Fig 1A and 1B).
In the susceptible C57BL/6J background, initial prenatal exposure to DEHP induced an intergenerational inheritance of anti-androgenic and/or pro-estrogenic symptoms, as well as a transgenerational inheritance of decreased sperm velocity affecting all generations tested up to F3 (Fig 1C and 1D). Importantly, the multigenerational impact found in the C57BL/6J exposed lineage cannot be formally proved in the present study. Indeed, the experimental design lacks F2 and F3 control conditions in both backgrounds, i.e. "FVB.CTL.F2", "C57.CTL.F2" and "C57.CTL.F3". This was due to the priority that was given to perform comparisons between both strains. However, the DEHP-induced transgenerational transmission of decreased sperm velocities recorded in the exposed C57BL/6J seems to be consistent when compared with stable levels recorded for the same parameters in the exposed lineage of the FVB/N strain (Fig 1B  and 1D). Previous reports strongly suggest that DEHP may trigger pro-estrogenic and antiandrogenic effects in a multigenerational manner, at least in some strains of mice and rats. Transgenerational inheritance induced by prenatal exposure to DEHP was previously established for cryptorchidism and testicular germ cell in CD-1 outbred mice and Sprague Dawley rats [128,129], for altered estrous cyclicity and decreased folliculogenesis in female CD-1 mice prenatally exposed to DEHP or to environmentally relevant phthalate mixture [130][131][132][133]. Targets for DEHP-susceptibility (Fig 3B and 3C) are the androgen-dependent Saa2, the serine peptidase inhibitor Kazal type 1 Spink1, the poorly described 9530003J23Rik, as well as the androgen-dependent Sva, Svs2, Svs3a, Svs3b, Svs4, Svs5, Svs6 and Pate4 genes (Table 1). Previous independent reports showed a downregulation of Svs5 in the fetal testes of Sprague-Dawley outbred CD rats exposed in utero to 750 mg/kg per day of DEHP during E12-19 in association with reduced AGD [134]. Targets were silenced by methylation of their promoters (Fig 3B) that contained AR binding sites, notably in Sva, Svs2 and Svs3ab (Fig 3C). In addition, the targets exerted a direct impact on sperm motility ( Table 1). The analysis revealed SNPs (rs28279710, rs46648903, rs46677594, rs32977910, rs48287999) affecting sexual hormonebinding sites specific to C57BL/6J and co-localized in the DNA encoding Svs2, Svs3ab, Svs4, Svs6 and Svs5 genes. The SNPs may explain the strain-specific DEHP response by rendering the allele more responsive to SHP in C57BL/6J compared with FVB/N (Fig 4C). In C57BL/6J, an apparent DEHP-induced transgenerational inheritance of decreased sperm velocity ( Fig  1D) was induced by in utero exposure to DEHP and correlates with both silencing of the Svs3ab gene (Fig 5B) and the decreased production of SEMG2 (Fig 5D). Among genes encoded by the susceptible allele, Svs2 encoding SEMG1 and Svs3ab encoding SEMG2 are responsible for the production of the two most abundant proteins in sperm known to regulate sperm motility [135].
In conclusion, prenatal exposure to the anti-androgenic DEHP molecule induces intergenerational and transgenerational alterations observed in a susceptible mice background, but not in a resistant mice strain. We found converging evidence for strain-specific, SNP-mediated alterations in androgen-regulated sperm transcripts in association with the alteration of male fertility parameters transmitted across several generations in C57BL/6J. This transmission may be mediated by the epigenetics of the paternal DNA having escaped the erasure-reprogramming system or by other spermatozoa-dependent unknown mechanisms. Biological support for transmission across generations may also involve dysregulated sperm RNA. However, injections of sperm RNA into fertilized oocytes, technically difficult, would be required to formally prove their transmission of the phenotype. We observed that strain-specific genomic variations may confer complex resistance mechanisms. SNPs were identified as associated with the strongest DEHP-induced specific changes in the germ cell content of RNAs coding for proteins required for sperm motility. Finally, the relevance of these findings for humans remains to be investigated. By extrapolation, the results suggest that human embryos exposed to DEHP may not be equally affected in their future reproductive health. According to our findings in mouse, unequal human susceptibility to prenatal exposure to DEHP should first imply genomic variations affecting directly SHP and SNP located in androgen responsive element (ARE). SNP located in ARE, as well as in FOXA1 and ERG binding sites in Human were enriched in prostate cancer susceptibility loci [136], whereas MEHP may advance the progression of prostate cancer through activating the hedgehog pathway [137]. If the purpose is to assess the impact of DEHP on Semenogelin production in human semen, SNP characterized as a quantitative trait loci (eQTL) for SEMG1 and/or SEMG2 should be taken into account. No eQTL are reported for SEMG2 according to the GTEx portal, whereas 181 eQTL were present in SEMG1. Among these eQTL, SNP Rs2233896 introduced a missense in SEMG1 affecting its expression with a minor allele frequency at 10%, and could therefore represent a good starting candidate. This study shows how mouse genomic variations may be associated with complex resistance mechanisms to prevent the inter-and transgenerational inheritance of altered male fertility parameters. More fundamentally, this study presents a good model to better characterize the array of endocrine signaling and target genes involved in the male reproductive function.

Materials and methods
This study was approved by the Ethics Committee for Animal Experimentation of the University of Geneva Medical School (Geneva, Switzerland) and by the Geneva Cantonal Veterinarian Office (permit reference: G61/3918) under the licenses GE/9/15 from February 2015 to August 2016 and GE/115/16 since 2017.
All resources used in the present work are given with the Resource Research Identifiers (RRID) number when possible, as recommended by NIH guidelines (https://scicrunch.org/). All the metadata generated in the study were deposited in a publically-available registry for transparency purposes. The data have also been deposited in the NCBI Gene Expression Omnibus [138] as GEO series accession number GSE107839. Files were converted by NCBI to Sequence Read Archive (SRA) files. We reused 10 samples from the previous GEO series accession number GSE86837. The MBD-seq data were reused from the GEO series accession number GSE67159. The analysis of SNPs impacting on hormonal binding motifs in the function of strain-specific RNA responses in male germ cells was deposited in the Mendeley repository (https://data.mendeley.com/datasets/3s94xbbtjx/draft?a=10be535a-1801-4805-977a-1b8d83b058f7). . DEHP is never mentioned in the compositions of these products and that material is not suspected to contain DEHP or other phthalates; DEHP is present principally in flexible PVC. However, we cannot exclude traces of DEHP or other phthalates in the animal Core Facility. If DEHP is present in the cages and in the bottles, its level should be very low, otherwise the material would be flexible. Moreover, the levels of DEHP contaminations through the equipment used at the animal Core Facility should be the same across the different experimental conditions, and cannot explain the phenotypical differences observed across these conditions. The ARRIVE checklist for this study is provided (S3 Table).

Prenatal exposure to DEHP and filial generations
One male and one female were mated in the same cage during the night and the presence of a copulatory plug the next day was considered as embryonic day 1 (E1). Prenatal exposure was performed daily from E9-19 covering the sexual differentiation period of the embryos by per os injections, using a micro-pipette, of 20 μl fixed volume into the mouth cavity of restrained pregnant female mice. Corn oil delivery vehicle for fat-soluble compounds (Sigma, C8267) was injected in controls (CTL) and 1.15 M DEHP (Fluka-Sigma, DEHP, 80030) diluted in corn oil was injected in exposed mice (D300). The dose of DEHP (CAS No. 117-81-7) was calculated for an estimated mouse weight of 30 g and corresponded to 300 mg of DEHP per kg of mice per day. F1 males prenatally exposed to DEHP were crossed with newly-ordered females of the respective backgrounds to produce F2. The third filial generation (F3) was obtained in the C57BL/6J background only by crossing F2 males with newlyordered females. The experimental design is shown (Fig 1A and 1C).

Computer-assisted sperm analysis and RNA extraction from spermatozoa
Sperm samples were extracted from 100 days of life age-standardized males. The left and right cauda epididymides of each mouse were cleanly dissected and deposited in a petri dish (Falcon, 351008) filled with a pre-warmed 1 ml M2 medium (Sigma, M7167). Epididymides were incised to allow the sperm to swim out during 10 min at 37˚C. Sperm suspensions were transferred in a 100 μm-deep homemade chamber placed on a stage warmer to maintain the temperature at 37˚C during the CASA process (Minitherm, Hamilton Thorme Lac, Beverly, MA, USA). Observations were made with a 4X phase contrast Objective (Olympus, RMS4X) at a final magnification of 40x with fixed parameters. For each sample, a minimum of 400 sperm cells tracks were captured at a rate of 25 images per sec. In parallel, 900 μl of sperm suspension were collected without tissue debris and centrifuge to pellet the sperm. RNA extraction was then performed using TRIZOL according to the manufacturer's recommendations and including glycogen as carrier (Invitrogen, UltraPure Glycogen, Carlsbad, CA, USA). Elution was in 20 μl water (Bioconcept, Water for Molecular Biology DNA/DNAse/RNAse free) and samples were snap-frozen and conserved at -80˚C.

Estimation of the purity of spermatozoa samples
The purity of spermatozoa used to extract RNA was assessed before the procedure by visualizing sperm under a contrast-phase microscope during the CASA process. After RNA sequencing, the purity of sperm RNAs was assessed by read length distributions as corresponding to the profile of mature sperm RNA signatures, involving a peak at 22 for microRNA and a peak at 32 nucleotide for tRNA-derived small RNAs [139]. The read length distribution of all samples is reported (S1 Fig).

Preparation of all-RNA-seq libraries and sequencing
All-RNA-seq was performed as previously described [13]. RNA molecules were treated to be compatible with adapter ligation following a protocol developed by Fasteris SA (Geneva, Switzerland), resulting in RNA molecules carrying a 5'-monophosphate extremity and a free 3'hydroxyl group independent of their starting differences (Fig 2A). Libraries were prepared using the TruSeq small RNA kit (Illumina Inc., San Diego, CA, USA) and sequenced on a HiSeq 2500 (Illumina Inc.). Base calling was performed using HiSeq Control Software 2.2.58, RTA 1.18.64.0 and CASAVA-1.8.2.

All-RNA-seq reads mapped to mm10 genome and differential quantification
Fastq files were processed on a Linux workstation (Ubuntu 14.04 LTS) and on the "Baobab" high performance computing cluster of the University of Geneva. The analysis was performed using TopHat and Cuffdiff [140]. Quality control of the sequencing process was performed with fastqc ensuring >Q28 in the quality control report for all reads. RA3 adapter sequence was clipped from fastq files and the minimum reads size was settled at 18 pb. Reads were mapped to the mm10 mouse genome using TopHat with appropriate strand-specificity options activated ($ tophat -p 16 -G genes.gtf-library-type fr-secondstrand -o genome reads. fastq). BAM files were indexed and sorted with samtools and reads statistics were computed to summarize the mapping process (S1 Table). Quantifications were performed using TopHat.

Hierarchical cluster analysis and heat map representation
Hierarchical cluster analysis was performed in R with iterative tests involving the diffData function of CummeRbund [141]. A two-dimensional graphical false-color image representation (heatmap) of the gene expression data was produced in R using the heatmap.2 function of the gplots package (Fig 2B) [142].

Analysis of pathways
Selected genes were tested for enrichments using STRING [143,144]. The analysis was performed online (http://string-db.org/) with the list of official gene symbols and selecting the "Mus musculus" organism. The minimum required interaction score was set at the highest confidence level (0.9). A significant enrichment was defined based on a false discovery rate (FDR) <0.01. The representative enrichment per category are reported (S2 Table).

Expression of identified targets in germ cell subtypes
Expression of genes across male germ cell subtypes came from a previously published independent study [15]. Primary spermatocytes and haploid round spermatids were isolated from B6129SF2/J mice at 8-9 weeks, whereas spermatogonia were isolated from testis of 4-8 days postpartum mice using enzymatic digestions of decapsulated testis combined with sedimentation techniques. In our previous study, enriched germ cell populations were analyzed using the microarrays "A-AFFY-45" (Affymetrix GeneChip Mouse Genome 430 2.0). The data were deposited in the ArrayExpress database archived in the European Bioinformatics Institute under the accession number "E-TABM-130". Data were uploaded and integrated in our analysis (Fig 2G).

Combined promoter methylation and RNA expression analysis
All-RNA-seq and MBD-seq data available for all four lineages at F1 were merged by gene symbol in R software separately for both strains (Fig 3B). For each gene, the promoter is defined as the region spanning 2000 bp upstream of the gene's 5'-end to 200 bp downstream. Its methylation level is measured as the number of reads obtained in a 2.2 kilobase probed region divided by the number of CpG in the tested region [11].

Detection of binding sites for the sexual hormone signaling performed in the promoters of the identified targets
Binding sites for FOXA1-3, ESR1 and AR were detected in promoters from -2000 to 0 bp among the selected targets using Rsat [111] (Fig 3B). Briefly, promoter sequences were extracted from Mus musculus GRCm38 using a list of official gene names and an analysis was performed using the matrix-scan full option with the same HOCOMOCO position weight matrices (PWM) used afterwards [145].

SNP-dependent hormonal binding motifs analyzed in function of strainspecific RNA responses to DEHP exposure
Bioinformatic approaches were used to assess strain-specific SNP affecting binding sites for hormonal controls and genes strain-specifically dysregulated by prenatal exposure to DEHP on genome-wide scale. The R packages and libraries "VariantAnnotation" [146], "Biostrings" [147], "BSgenome.Mmusculus.UCSC.mm10" [148] and "EnsDb.Mmusculus.v79" [149] were required for the analysis. 5'556'605 FVB/NJ SNP variations compared with C57BL/6J were uploaded from the FTP site of the Sanger Institute and filtered to 4'823'906 of high quality. Filtrated SNPs between both strains were used to reconstruct C57BL/6J and FVB/NJ strain-specific alleles with sizes of sequences equal to two times the size of HOCOMOCO position weight matrices (PWM) for mouse AR (named ANDR in the HOCOMOCO PWM and replaced by AR in the present text for simplification purpose), ESR1 and FOXA1-3 [145]. The "matchPWM" function in R was used to computed strain-specific identification scores based on the SNP and the respective scores were reported for the reference (C57BL/6J = ref) and the alternative (FVB/NJ = alt) alleles in both the reverse complement and the positive strand [150]. The results were annotated reporting the distance and the name of the most proximal genes from the SNP-dependent motifs. Genome-wide statistical significances of SNP-dependent detected motifs were based on the HOCOMOCO threshold to p-value relationships. Tables were created with the resulting detected SNP-dependent motifs. Generated tables were then merged individually to the all-RNA-seq derived FPKM data and deposited in Mendeley. SNPs affecting the expression of genes in a strain-specific manner by abolishing or creating a binding motif for FOXA1-3, AR and ESR1 were identified using δ, for motif specificity, in function of Z, reflecting strain-specific transcriptional changes; (delta <-1 or  (Table 2).

Validation experiments involving bisulfite-pyrosequencing
DNA was extracted using ethanol precipitation from the organic phase of TRIZOL conserved at -20˚C. These samples correspond to the sperm samples obtained during the CASA process that were used for RNA extraction. Bisulfite conversion of DNA was performed using the EZ DNA methylation-Lightning kit (Zymo Research, Irvine, CA, USA; D5030). All assays were validated using methylation standards (EpigenDX, Hopkinton. MA, USA; 80-8060M-PreMix). PCR cycling conditions were 95˚C for 15 min, then 50 cycles of 94˚C for 30 s, 54˚C for 30 s (replaced by 56˚C for 30 s for Svs3b), 72˚C for 30 s, with a final elongation during 1 min at 72˚C. Pyrosequencing was performed using a Pyromark Q24 instrument and a workstation from Qiagen (Hilden, Germany) using appropriated enzyme, substrate, nucleotides and reagents. For methylation measured in the Svs2 promoter, a 175 bp amplicon encompassing 2 CpG sites was generated using forward 5'-AGGAGGGGTTAGTGTTTTTTGTAGT-3' and biotinylated reverse 5'-[BIO] CAACCACAACTAACTATCTTTCCAAA-3' primers. The sequencing primer was 5'-GTTATATTTAGTTAGAGAAGATATAAA. The sequence to analyze was YGTGTTAGTTAGTTTAGTTTTTAGTGAAGGTTTTTTTTGATAAGATGAAGTTTTTTGT TTTYGTTTTTTTTTTGTTTTTTATTTTGGAAAGATAGTTAGTTGTGGTTG and the dispensation order was GTCGTGTAGTAGTAGTTAGTGAGTTGACTAGATGAGTTAGTTCTGTTG. The two CpG sites were validated by methylation standards. For methylation measured in the Svs3b gene, a 118 bp amplicon encompassing 2 CpG sites was generated using the forward primer Svs3b_F 5'-GGGAGAATGTTTATAAGGATGTTATG-3' and the reverse Svs3b_R_biot 5'-[BIO]CATCAACCCACTACTATACCCAAAC3'. The sequencing primer was the forward primer. The sequence to analyze was GTAGTGAGAGAATTTTTTAYGGGGTGATTTTTTGGAT GGAYGTTGTTTTTTTTATTTATATTTATATGTTTGGGTATAGTAGTGGGTTGATG and the dispensation order was AGCTAGTGAGAGATTGATCGGTGATTGATGTATCGTG. Both CpG sites were validated using methylation standards.

Western blots of SEMG1 and SEMG2
Sperm total protein samples were obtained from TRIZOL fractions after RNA and DNA extractions according to the manufacturer's recommendations by precipitation with 100% isopropanol.

Statistical analysis
One-way analysis of variance with post-hoc Tukey honestly significant difference test (ANO-VA-Tukey HSD) was used to identify means that were significantly different between conditions for the parameters measured with CASA, as well as for AGD and testes weight measurements (Fig 1). The data obtained during the CASA process were tested for following a normal distribution using the Shapiro-Wilk normality test, as well as for homogeneity of the variances across the seven experimental conditions using the Bartlett test of homogeneity of variances. The significance level was set at p � 0.05. In all-RNA-seq analysis, the CuffDiff computed false-discovery rate adjusted p-values (q-values) were used to identify the transcript levels statistically significantly different between the labeled conditions. Probability estimates for each cluster resulting from the hierarchical cluster analysis were obtained using Pvclust [151]. Pearson's chi-square goodness-of-fit test was computed in R version 3.4.2 using the "chisq. test" function to assess the statistical significance of differences in frequencies of the respective proportions of genes whose expression increases during spermatogenesis among the identified targets versus all the genes. Correlations between MBD-seq and RNA-seq data were analyzed using non-parametric Kendall's tau and Spearman's rho rank correlations in R (Fig 3B). Nonparametric pairwise Wilcoxon test with Benjamini & Hochberg correction for multiple testing was used to estimate differences of transcript levels between conditions associated with SNPs (Figs 3A and 4) [152]. Parametric Pearson correlations between CpG methylation and the ordered conditions involving C57BL/6J strain were performed in R (Fig 5B).
Supporting information S1 Fig. Read lengths distribution across samples. A Peaks specific for sperm RNAs were identified. The peak at 22 nucleotides is specific to the size of mature microRNAs. The peak at 32 nucleotides is specific to the size of tRNA-derived small RNAs. The peak at 50 bp involves mainly coding RNAs. (TIF)

S2 Fig. Example of SNP validation with all RNA-seq reads in a region of interest.
Rs28279704 missense mutation validated in the Svs4 gene as being "A" in C57BL/6J, reversely in the ATG codon for methionine, and "G" in FVB/N, reversely in the CG codon for threonine.