Assessing the Ability of Developmentally Precocious Estrogen Signaling to Recapitulate Ovarian Transcriptomes and Follicle Dynamics in Alligators from a Contaminated Lake

Background: Concern has grown in recent decades over anthropogenic contaminants that interfere with the functioning of endocrine hormones. However, mechanisms connecting developmental processes to pathologies associated with endocrine-disrupting chemical (EDC) exposure are poorly understood in naturally exposed populations. Objectives: We sought to a) characterize divergence in ovarian transcriptomic and follicular profiles between alligators originating from a historically EDC-contaminated site, Lake Apopka, and a reference site; b) test the ability of developmentally precocious estrogen exposure to recapitulate site-associated patterns of divergence; and c) test whether treatment with exogenous follicle-stimulating hormone (FSH) is capable of rescuing phenotypes associated with contaminant exposure and/or embryonic estrogen treatment. Methods: Alligators eggs were collected from a contaminated site and a reference site, and a subset of eggs from the reference site were treated with estradiol (E2) during embryonic development prior to gonadal differentiation. After hatching, alligators were raised under controlled laboratory settings for 5 months. Juveniles from both sites were divided and treated with exogenous FSH. Histological analyses and RNA-sequencing were conducted to characterize divergence in ovarian follicle dynamics and transcriptomes between sites, between reference and E2-treated animals, and between FSH-treated and nontreated animals. Results: We observed broad site-of-origin divergence in ovarian transcriptomes and reductions in ovarian follicle density between juvenile alligators from Lake Apopka and the reference site. Treating embryos from the reference site with E2 overwhelmingly recapitulated transcriptional and histological alterations observed in Lake Apopka juveniles. Ovarian phenotypes observed in Lake Apopka alligators or resulting from estrogen treatment were only partially rescued by treatment with exogenous FSH. Discussion: Recapitulation of ovarian abnormalities by precocious E2 revealed a relatively simple mechanism underlying contaminant-induced pathologies in a historical example of environmental endocrine disruption. Findings reported here support a model where the developmental timing of estrogen signaling has the potential to permanently alter ovarian organization and function. https://doi.org/10.1289/EHP6627


Introduction
The last century of human industrial development has been characterized by widespread production and release of anthropogenic chemical pollutants into natural systems. In the last three decades, mounting concern has developed regarding the ability of a subset of these environmental pollutants to interfere with the normal functioning of endocrine hormones and which are termed endocrinedisrupting chemicals (EDCs). The mechanisms of action underlying the biological impacts of exposure to EDCs range from hormone mimicry, receptor antagonism, alterations to hormone synthesis or clearance, and altered hormone transport throughout the body (Gore et al. 2015). Of these, the proclivity of contaminants to disrupt estrogen signaling is among the best studied modes of action for EDCs (Guillette 2006). Because estrogens regulate patterns of proliferation and differentiation during embryonic stages (Bondesson et al. 2015), inappropriate activation of signaling can induce dramatic shifts in tissue patterning that persist into later life. Consistent with Barker's developmental origins of health and disease model (Barker 2004), EDC exposure during development has been linked to numerous reproductive disorders in both women (Crain et al. 2008;Gore et al. 2015) and men (Gore et al. 2015;Sweeney et al. 2015).
Despite the fundamental utility of laboratory models for investigating mechanisms behind EDC-induced pathologies, these reductionist approaches fail to adequately capture environmental complexity. Consistent with the sentinel species concept (National Research Council 1991), wildlife species exposed in situ have repeatedly served as models for environmental health to fill this gap (Blus et al. 1974;Hickey and Anderson 1968;Gilbertson and Fox 1977;Guillette et al. 1994). Among wildlife models, studies on a population of American alligators inhabiting a highly contaminated freshwater system, Lake Apopka (AP), in central Florida (Orange County, FL) helped to pioneer our early understanding of how contaminant exposure during embryonic stages determines future reproductive health. There, developmental exposures to a complex mixture of organochlorine pesticides (OCPs) (Heinz et al. 1991;Rauschenberger et al. 2007Rauschenberger et al. , 2009) have been associated with a suite of reproductive abnormalities, including ovarian and testicular morphological abnormalities (Guillette et al. 1994) and shifts in gonadal transcription that persist into juvenile stages Kohno et al. 2008;Moore et al. 2012). Reproductive pathologies at AP have been repeatedly observed in juvenile animals collected as eggs and raised under controlled laboratory settings, demonstrating a persistent site-of-origin effect that appears to have developmental origins (Guillette et al. 1994;Hale et al. 2019;Moore et al. 2012). However, previous studies in this system, like many other ecological and epidemiological approaches, have relied on correlative methods to associate exposure to phenotype. Many OCPs that have been shown to specifically contaminate AP are capable of inducing estrogen receptor (ER) activation when assayed in vitro or to alter sex determination (Klotz et al. 1996;Matter et al. 1998;Milnes et al. 2005b;Vonier et al. 1996), leading us to hypothesize that the maternal deposition of these contaminants in alligator egg yolk might constitute an aberrant estrogenic signal during early gonadal development that ultimately contributes to altered ovarian function later in life.
In the present study, we leverage the ecological relevance of the alligator to investigate the functional consequences of developmental EDC exposures in naturally exposed organisms and to probe the causal mechanisms underlying complex contaminantinduced pathologies. Alligators are oviparous, and field-collected eggs can yield laboratory-reared offspring that are exposed to maternally deposited contaminants in ovo but which subsequently do not experience direct exposure through their postnatal environment, thus permitting investigations into the functional consequences of embryonic exposures that persist into later life stages.

Methods
Animal Husbandry, Follicle-Stimulating Hormone Administration, and Sample Collection All experiments described herein were reviewed and approved by the institutional animal care and use committee at the Medical University of South Carolina (Charleston, SC), and alligator egg collections were approved and permitted by the Florida Fish and Wildlife Conservation Commission. The husbandry for animals used in this study has been previously described in detail Hale et al. 2019). Briefly, alligator eggs were collected shortly after oviposition from two sites in Florida, AP (Orange county, FL) and Lake Woodruff National Wildlife Refuge (WO; Volusia county, FL). Both locations and populations are well characterized; AP is characterized by high levels of OCP contaminants, the product of long-term agricultural runoff and an acute spill event of dicofol, a dichlorodiphenyltrichloroethane (DDT) and DDTmetabolite-containing miticide, in the 1980s. These contaminants have been detected at high levels in adult and juvenile alligators and in egg yolk (Guillette et al. 1994(Guillette et al. , 1999Heinz et al. 1991;Rauschenberger et al. 2007Rauschenberger et al. , 2009Woodward et al. 2011). In contrast, WO is a site for which minimal anthropogenic disturbance or contaminant input has been detected and is used herein as a reference site (Crain et al. 1998;Guillette et al. 1994;Milnes et al. 2005aMilnes et al. , 2008Moore et al. 2010bMoore et al. , 2010cPickford et al. 2000). Eggs were collected from both sites in June 2014, candled to assess viability, staged according to Ferguson (1985), and transferred to artificial nests of damp sphagnum moss at the Hollings Marine Laboratory (Charleston, SC). For this study, 18 eggs total were collected from 5 clutches at AP; at WO, 31 eggs total were collected from 12 clutches. All eggs were collected at or prior to Stage 16. Eggs were maintained at 32°C until Stage 19, at which time, WO eggs were randomly distributed among two treatment groups and dosed once with either 0:5 lg=g (egg weight) estradiol-17b (E 2 ; Sigma-Aldrich) or vehicle (VEH) control (95% ethanol; Sigma-Aldrich) via a single topical administration to the apical surface of the eggshell. All AP eggs in this study were treated with VEH only (i.e., no AP-E 2 -treated embryos were used). Stage 19 of alligator development precedes the onset of aromatase expression , thus, exposure to exogenous estrogen or estrogenic compounds at this stage is referred to as precocious given that it precedes endogenous estrogen biosynthesis.
Following treatment, embryos were incubated at 30°C, a temperature that produces exclusively females. Upon hatching, the animals were individually marked with numbered monel tags between the middle digits on both hindlimbs and transferred to indoor fiberglass tanks at the Hollings Marine Laboratory that permit both basking and swimming. Neonates were maintained on 12 h/12 h light:dark cycles and fed ad libitum with a commercially available crocodilian chow (Mazuri Exotic Animal Nutrition) for approximately 2 months following hatching. To ensure that neonates were housed with similarly sized individuals, all animals were weighed and sorted by body mass every 2 wk. At 2 months of age, the animals were switched to a feeding schedule based on body mass to ensure relative homogeneity of animal size . Briefly, animals weighing <92 g were fed daily, animals weighing 92-132 g were fed three times per week, and animals weighing >132 g were fed twice weekly. At 5 months of age, animals across all treatment groups were further divided into two treatment groups and were dosed once daily for 4 d with either 277 lU=g recombinant ovine follicle-stimulating hormone (FSH; Sigma-Aldrich F8174) or VEH control (0.8% sterile saline) via intramuscular injections at the base of the tail. On the fifth day following initial treatment, animals were assigned sequential dissection identifiers and euthanized with 0:1 mg=g pentobarbital (Sigma-Aldrich) followed by decapitation. Ovaries were then necropsied and weighed; the right ovary was fixed in RNAlater ® , rocked for 12 h on an orbital shaker at 4°C, then frozen at −80 C. The left ovary was fixed in 4% formaldehyde (10% neutral buffered formalin) for 24 h and stored at 4°C in 70% ethanol. Total RNA was isolated from RNAlater ® -fixed right ovaries using a modified acid guanidinium thiocyanate-phenol-chloroform (AGPC) extraction with silica column purification for RNA sequencing (RNAseq) analysis. Briefly, gonadal tissues ( ∼ 10 mg per sample) were manually homogenized in an AGPC lysis solution (water-saturated acidic phenol, 2 M guanidinium, thiocyanate, 95 mM sodium acetate, 12 mM sodium citrate, 0.24% N-lauroyl sarcosine, 14:4 M beta-mercaptoethanol) with a Retsch ball mill and stainless steel beads at 30 Hz for 2 min. Following phase separation with 0:2 mL 37% chloroform, RNA-containing aqueous layers were mixed 1:1 with 100% ethanol and applied to a silica-membrane spin column (EconoSpin; Epoch Life Science). On-column DNase I digestion was conducted to eliminate genomic DNA contamination (5Prime DNase I), followed by elution in ultra-pure diethylpyrocarbonatetreated water. Final group sample sizes are reported in Table S1. All downstream analyses utilized identifiers assigned at dissection, thus authors were blinded as to treatment group or site of origin.

RNAseq Read and Count Data Generation
Alligator RNAseq library preparation and sequencing was conducted by the Georgia Genomics and Bioinformatics Core at the University of Georgia. Briefly, 1-2 lg total RNA (n = 7 libraries per treatment group, 42 libraries total) were assessed for concentration and quality via Agilent 2100 Bioanalyzer [RNA integrity number (RIN) ≥7:80, average = 8:6) and then diluted to a common concentration. Samples chosen for library preparation were selected randomly. Libraries were prepared from poly(A)enriched mRNA using the KAPA Stranded mRNA-Seq kit for Illumina platforms (KAPA Biosystems) and were sequenced on an Illumina NextSeq (150 cycles, 75-bp paired-end reads) in two individual sequencing runs. Sample libraries used to assess population-level transcriptome divergence and FSH responsiveness were prepared and sequenced simultaneously (i.e., WO-FSH, WO-VEH, AP-FSH, and AP-VEH). Sample libraries used in assessment of the effects of precocious estrogen signaling were prepared and sequenced independently of the previous run (i.e., WO-VEH and WO-E 2 -VEH). Identical WO-VEH samples were sequenced in both experiments, but libraries were prepared independently (i.e., libraries were not resequenced but, rather, were generated de novo). Read quality was assessed via FastQC (https://www.bioinformatics.babraham.ac.uk/projects/ fastqc/) and MultiQC (version 1.5) (Ewels et al. 2016). Alignment summaries and library sizes for each experiment are reported in the supplemental material in "Population-level transcriptome divergence and FSH responsiveness" and "Effects of precocious estrogen signaling." Following quality assessment, reads from both sequencing runs were determined to lack adapter contamination (<0:1% in all samples) and exhibited mean read quality Phred scores >30 and were therefore not subjected to quality trimming. Fastq reads were aligned to the most recent alligator genome assembly (ASM28112v4) (Rice et al. 2017) using Hisat2 (version 2.1.0) (Kim et al. 2015) with parameter -dta and guidance parameters -ss and -exon enabled during index generation; splice site and exon coordinate information were generated from alligator National Center for Biotechnology Information Reference Sequence (NCBI Refseq) annotations (GCF_000281125.3, annotation release 102) using BEDtools (Quinlan and Hall 2010). The resulting alignments were coordinate sorted and converted to .bam format via Samtools (version 1.6) (Li et al. 2009). Once sorted, count matrices were generated using packages GenomicFeatures and GenomicAlignments (Lawrence et al. 2013) in R (version 3.5.1; R Development Core Team). The GenomicFeatures function makeTxDbFromGFF was used to generate exon-by-gene coordinates for the alligator assembly, which were then used in GenomicAlignments to generate feature counts from sorted alignments, using the function summarizeOverlaps (parameters mode = Union, fragments = TRUE, singleEnd = FALSE, ignore:strand = FALSE). This generated a count matrix with 24,848 unique genes in both sequencing runs.

Histology and Ovarian Morphology
Formaldehyde-fixed ovaries were bisected symmetrically on their transverse midline so as to evenly split the ovaries into two equally sized halves, which were paraffin embedded and sectioned at 4 lm thickness. Sections were hematoxylin and eosin stained and imaged, then composite-stitched using a Keyence BZ-X710 and Keyence automated XYZ image stitching software at 10 × and 20 × magnification. Images were analyzed in ImageJ (version 1.52a; Schneider et al. 2012). Margins of the ovarian cortex was identified by its distinct border with underlying ovarian lacunae and trabeculae, and these landmarks were used to ensure comparable ovarian regions were compared across samples. Total cortical area was measured using the freehand selection tool. Stage III oocytes were identified by a pink basophilic cytoplasm and by the presence of a complete follicular granulosa cell layer and one or more surrounding thecal cells (Moore et al. 2010a;Uribe and Guillette 2000). Follicle diameter was measured for each oocyte on its longest axis. When possible, images from both ovarian halves were analyzed, and the resulting measures averaged to represent each ovary; one representative section or pair of sections were analyzed from each sample. Follicle density is reported as the number of late Stage II (oocytes with pink basophilic cytoplasm that lacked a definitive follicular granulosa layer and thecal cells) and Stage III oocytes normalized to total cortical area.

Differential Gene Expression and Follicle Analysis
Identification of differentially expressed genes was conducted using package edgeR (Chen et al. 2016;McCarthy et al. 2012;Robinson et al. 2010;Robinson and Smyth 2007) in R (version 3.5.1; R Development Core Team). In the first experiment describing population-level effects between AP and WO (no E 2 -treated animals), low-expression genes [count per million (CPM) <1, minimum number of expressing libraries = 7] were removed prior to analysis (with 18,435 genes passing filtering). In the second experiment comparing WO controls and WO-E 2 animals, library sizes were approximately doubled, thus filtering was relaxed (CPM <0:5, minimum number of expressing libraries = 7), retaining 19,230 genes for analysis. Multidimensional scaling and principal components analysis was used to identify and remove two WO-FSH-treated samples collected from two different clutches because of their relative distance from remaining libraries in that treatment group and high biological coefficient of variation (BCV) values ( Figure S1; n = 5 after removal; Samples 60 and 68). Removal of these two outlier samples reduced the common dispersion value from 0.09 (BCV = 0:30) to 0.059 (BCV = 0:24), which concomitantly reduced a band of relatively high BCV genes observed by plotting log(CPM) against BCV ( Figure S1B,F). Histological examination of these samples revealed that they exhibited exceedingly low Stage III follicle counts in comparison with other WO samples; thus, they were excluded from transcriptomic and histological (described in following paragraph) analyses. Following filtering in both experiments, libraries were trimmed mean of M-values (TMM) normalized to adjust for composition biases and fit to a negative binomial model (glmQLFit, robust = TRUE). Hypothesis testing was conducted with planned linear contrasts via quasi-likelihood F-tests (glmQLFTest). Genes with a Benjamini-Hochberg adjusted false discovery rate (FDR; Benjamini and Hochberg 1995) p < 0:05 were considered significant. Overlap of differentially expressed genes (DEGs) between WO and AP basal contrasts and WO-E 2 -exposure was assessed by identifying presence/absence of significantly affected genes across experiments; significant recapitulation of genes elevated and suppressed in AP and in WO E 2 -exposed animals was determined using the hypergeometric distribution (a = 0:05). Because total features passing low-expression filters differed slightly between experiments (n = 18,435 vs. n = 19,230), only shared features passing filtering in both experiments were used in the recapitulation analyses (n = 18,360 shared genes). Overlap between gene lists was identified using Venny (http://bioinfogp.cnb.csic.es/tools/venny/index.html).
Total body mass at necropsy (in grams) was analyzed via two-way analysis of variance (ANOVA), and follicle density (counts per unit cortical area) data were analyzed with one-way ANOVA in GraphPad Prism (version 8.0.1). Homogeneity of variances across treatment groups was confirmed via Spearman (two-way ANOVA) or Brown-Forsythe (one-way ANOVA) tests (a = 0:05). Residual normality was confirmed via the Shapiro-Wilk test of normality (a = 0:05). All treatment groups for the two metrics assessed met both assumptions. Post hoc tests were conducted via Tukey's multiple comparisons test within FSH groups (FSH-treated and nonchallenged VEH groups, independently). Rescue of ovarian transcriptomes by FSH, defined as a reduction in number of DEGs between populations, was assessed via Fisher's exact test (a = 0:05), comparing the relative proportions of genes that were not differentially expressed between FSH-treated groups but that were differentially expressed in initial contrasts (non-FSH-treated groups, AP-VEH vs. WO-VEH) (i.e., a significant reduction in the proportion of DEGs is indicative of rescue).
RNAseq read counts and qPCR expression values were assessed individually using two-way ANOVA (differences in expression between FSH-treated and VEH-treated juveniles from both AP and WO) or unpaired t-tests (basal differences in expression between WO-VEH-treated and AP-VEH-treated animals) (a = 0:05). Significance of effects were then compared across quantification approaches. Residual normality and homoscedasticity were confirmed via the Shapiro-Wilk test (a = 0:05) and Spearman's test (a = 0:05), respectively, for both approaches. Values were transformed as necessary to meet these model assumptions.

Annotation of Uncharacterized Loci
During DEG analyses, we observed a large proportion of uncharacterized loci in RefSeq annotations (5,469 genes of 18,435 total passing filtering in the first sequencing experiment and 6,077 of 19,230 in the second), many of which were consequential for functional enrichment [e.g., CYP19A1 (LOC102566432), CYP17A1 (LOC102567971), and CYP11A1 (LOC102569028)]. To confirm the identity of these genes, read alignments were assembled into transcripts using StringTie (version 1.3.3) in a two-pass assembly approach (Pertea et al. 2015(Pertea et al. , 2016, merging transcripts across all libraries (parameter -merge). Fasta sequences were then extracted from the StringTie-merged assembly via the GFFread utility in Cufflinks (version 2.2.1) (Trapnell et al. 2010) and used in NCBI Blastx (https://blast.ncbi.nlm.nih.gov/ Blast.cgi) against the Uniprot Swiss-protein database (UniProt Consortium 2017) (e-value cutoff = 1 × 10 −5 ). Top hits according to e-value for each unannotated transcript were used as evidence for gene identity in downstream analyses. This approach recovered identities for 2,666 genes (48.7% of the unidentified genes) in the population divergence and FSH response experiment, and 2,757 (45.4% of the unidentified genes) in the precocious E 2 experiment. Any LOC loci that was not identified via this approach was excluded from downstream analyses.

ER Response Element Enrichment of DEGs
The potential for direct involvement of estrogen signaling in driving transcriptome divergence was assessed by identifying genes in relative proximity to an estrogen response element (ERE). Briefly, alligator genomic regions matching the full palindromic human ESR1-binding motif (Lin et al. 2007) were predicted in silico using PoSSuM-search (Beckstette et al. 2006) (parameters -esa, -lazy enabled) using a p cutoff of 4:39 × 10 −6 , as described by Rice et al. (2017). This generated a list of 5,500 putative EREs in the alligator genome. In order to capture both proximal and highly distal regulatory interactions between predicted EREs and genes, three different distance cutoffs separating a putative ERE and a coding sequence were used to define ERE-proximate genes: 250, 50, and 5 kb. Fisher's exact tests were then used to assess whether the proportion of ERE-proximate genes at each distance cutoff in DEG lists differed relative to the proportion of ERE-proximate genes in a background consisting of all detectable ovarian genes.

Functional Enrichment of DEGs: Transcription Factors, Gene Ontology, and Kyoto Encyclopedia of Genes and Genomes Enrichment
Functional annotation of Gene Ontology (GO) terms biological process, Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways, and enriched transcription factors for affected genes was conducted using g:Profiler (Reimand et al. 2016) and Enrichr (Chen et al. 2013;Kuleshov et al. 2016). Analyses were limited to genes changing by at least 2-fold between treatment groups. Transcription-factor enrichment was assessed using the ChIP-X Enrichment Analysis (ChEA 2016) database (Lachmann et al. 2010) through Enrichr, using an adjusted p < 0:05 to identify significantly enriched terms. To confirm these results in an alligator ovary-specific context, enrichment in highly affected genes lists [DEGs at least 2-fold responsive; log 2 fold change ðFCÞ > j1j] was compared with total enrichment in all annotated ovarian genes passing filtering (background) using Fisher's exact tests (a = 0:05) [ Table S3 (AP-responsive) and Table S5 (E 2 -responsive)]. GO and KEGG pathway analysis was conducted using unordered enrichment testing in g:Profiler (version r1760_e93_eg40), with a Benjamini-Hochberg adjusted FDR p-value correction and strong hierarchical filtering against an alligator ovary-specific background [with all genes passing filtering excluding nonannotated features (LOC genes with identities not retrieved through Blastx)]. Prior to enrichment testing, alligator annotations were converted to human Ensembl gene identifiers. Electronic annotations were excluded.

Gene Expression Clustering and Trait Association Analyses
Associations between gene expression patterns and follicle densities were assessed through pairwise Pearson correlations using the package psych (version 1.8.12) (Revelle 2018) in R (version 3.5.1; R Development Core Team) in RStudio. Read counts passing filtering in both experiments were variance-stabilizing transformation (vst)-transformed using the package DESeq2 (Love et al. 2018). Pearson correlation coefficients between sample follicle densities and expression values were then determined using the psych package function corr.test with a Benjamini-Hochberg FDR-adjusted p-value (a = 0:05).

ER Expression in the Embryonic Gonad
Ontogeny of ER expression in the developing alligator ovary was investigated using RNA-seq reads accessed through the NCBI Sequence Read Archive (SRA; accession: PRJNA322197). Briefly, reads were sourced from alligator embryonic gonads at Stages 19, 20 (Stage 19 + 3 d), and 26 (Stage 19 + 30 d) incubated at a femalepromoting temperature (Rice et al. 2017). Resulting reads were trimmed using TrimGalore (https://github.com/FelixKrueger/ TrimGalore) with the adapter overlap parameter -stringency set to 3 and with default Phred score quality and post-trimming read length parameters. Following trimming, validated paired reads and reads orphaned by trimming were each aligned independently to the alligator genome, as described above (alligator genome assembly ASM28112v4). Alignments were then sorted, .bam converted, and used to generate read counts, as described above. Read counts from both paired read and orphan alignments were summed to generate total counts for a given gene and then CPM transformed in edgeR.

Results
To address persistent shifts in ovarian function across populations, we investigated transcriptomic divergence in ovaries from Significantly enriched transcription factors predicted to regulate highly [log 2 ðFCÞ > j1j] suppressed genes at AP relative to WO controls. Significant enrichment for specific transcription factors was not detected for highly up-regulated genes at AP. Note: BMI1, BMI1 proto-oncogene, polycomb ring finger; BMP, bone morphogenetic protein; BP, biological process; CBX2, chromobox 2; ECM, extracellular matrix; ESR1, estrogen receptor 1; EZH2, enhancer of zeste homolog 2; FC, fold change; FDR, false discovery rate; GO, Gene Ontology; EED, embryonic ectoderm development; JARID2, Jumonji and A T-rich interaction domain containing 2; KEGG, Kyoto Encyclopedia of Genes and Genomes; MTF2, metal response element binding transcription factor 2; PHC1, polyhomeotic homolog 1; REST, RE1 silencing transcription factor; RING1B, ring finger protein 1b; RNAseq, RNA sequencing; RNF2, ring finger protein 2; SRP, signal-recognition particle; SUZ12, polycomb repressive complex 2 subunit; TP53, tumor protein p53; ZFP281, zinc finger protein 281. juvenile alligators originating from AP and WO, a commonly used reference site in this system because of its geographical proximity and minimal anthropogenic disturbance ( Figure 1A). We identified an unexpectedly large population effect ( Figure 1B,C), with 75.4% of detectable genes differing significantly between sites (FDR p < 0:05). Functional enrichment of strongly suppressed genes at AP (≥2-fold reduction relative to WO animals) revealed affected pathways related to cell cycle progression and proliferation ( Figure 1D; Excel Table S1) as well as pathways contributing to development and morphogenesis. In contrast, strongly upregulated genes at AP were enriched for ciliary and cytoskeletal terms as well as for extracellular matrix-related processes (Excel Table S2).
Given the large number of DEGs across sites, we investigated whether specific transcriptional regulators were predicted to disproportionately target affected genes. Our analysis revealed strong enrichment for polycomb group proteins ( Figure 1E; Table S3), particularly components of polycomb repressive complex 2 (PRC2). This included core components SUZ12, EZH2, and EED, as well as accessory factors JARID2 and MTF2. Consistent with the role of PRCs in establishing repressive transcriptional patterns, this enrichment was only observed in genes suppressed at AP. Interestingly, suppressed genes were also significantly enriched for putative targets of ERa (ERa=ESR1).

Recapitulation of AP Ovarian Transcriptomes by Developmentally Precocious Estrogen Signals
We next sought to explore potential causal mechanisms underlying the observed population-level divergence of the ovarian transcriptome. To test this, we exposed WO alligators to E 2 at Stage 19 ( Figure 2A), a stage that precedes the detection of aromatase expression and sexually dimorphic gonadal steroidogenesis Parrott et al. 2014). Using data from a study of alligator gonadal development (Rice et al. 2017, NCBI SRA accession: PRJNA322197) we observed that both estrogen receptors, ESR1 and ESR2, are expressed at this stage and at a level comparable to later stages of development ( Figure 2B; numerical figure summary in Table S4), suggesting that precocious estrogenic cues could directly induce canonical estrogen signaling in the early gonad. Thus, we directly assessed the ability of a precocious estrogenic cue to recapitulate AP transcriptional programming in the ovary. Similar to patterns observed at the population level, we observed a large degree of transcriptional divergence in WO-E 2 animals when compared with VEH-treated WO animals ( Figure  2C), with over 70% of genes differentially expressed between the two groups (FDR p < 0:05). Strikingly, estrogen treatment was capable of significantly recapitulating AP transcriptional programs in reference animals when compared with their control counterparts, inducing persistent up-regulation of over 76% of the same genes elevated at AP (relative to WO) and down-regulation of over 77% of genes with decreased expression in AP alligators (relative to WO) ( Figure 2D). These patterns remained significant in highly responsive (>2-fold affected) genes (39.3% highly up-regulated, 57.7% highly down-regulated; Figure S2A). Furthermore, targets for the same transcriptional regulators identified at AP (i.e., PRC components) were enriched across genes that were highly suppressed by E 2 treatment ( Figure S2C; Table S5). Similar to the ovarian transcriptomes of AP alligators, we also detected functional enrichment of cell cycle pathways in down-regulated genes as well as recapitulation of ciliary functional terms in elevated genes in ovaries from E 2 -treated WO animals ( Figure S2B).
Given evidence that OCPs were acting through an estrogenic mode of action, we sought to further explore the possibility of a role for the ER in mediating transcriptional divergence. To test this hypothesis, we identified differentially expressed genes [log 2 FC >j1j] in varying degrees of proximity to predicted EREs in the alligator genome. We found that genes with significantly greater expression in AP or E 2 -treated WO animals were consistently more likely to be proximate to an ERE when compared with all detectable ovarian genes. Interestingly, this pattern was reversed for genes with lower expression in AP and WO-E 2 animals as they were less likely than expected to lie in proximity to an ERE relative to background ( Figure 2E; numerical figure summary in Table S6).

Altered Ovarian Follicular Profiles
We next sought to investigate follicular profiles in AP and E 2 -treated WO animals (Figure 3). The majority of follicles observed were definitively Stage III (85.9%) and late Stage II were rare. Thus, both stages were combined into a single measure of follicle count per unit cortical area. Furthermore, average follicle diameter did not differ significantly by site ( Figure S3), indicating that the groups did not differ in relative proportion of these two follicle types, nor did total body mass at necropsy ( Figure S4). Interestingly, follicle density differed significantly between groups (p = 0:0015, F 2,19 = 9:371); the relative density of previtellogenic, Stage III follicles was markedly lower in AP animals compared with WO controls (AP vs. WO Tukey-adjusted p = 0:0011) ( Figure 3A-D). Developmental treatment of WO embryos with E 2 was sufficient to phenocopy the lowered follicle density observed in animals from AP given that both AP and E 2 -treated WO animals (WO-E 2 vs. WO Tukey-adjusted p = 0:0404) exhibited significantly lower density of Stage III follicles than WO controls but were indistinguishable from one another (AP vs. WO-E 2 Tukey-adjusted p = 0:4037) ( Figure 3E).
To better understand the relationship between altered follicular dynamics and the ovarian transcriptome, we employed pairwise correlations between follicle density and expression values for each gene detected, and we subsequently observed that the bulk of differentially expressed genes were highly correlated with follicle density (Figure 3F-G; numerical figure summaries in Table S7). Consistent with patterns of follicle density, the majority of genes with elevated expression [log 2 ðFCÞ < -1] in AP and WO-E 2 -treated animals compared with WO animals were negatively associated with follicle density, whereas most genes with reduced expression in these groups [log 2 ðFCÞ > 1] were positively associated with follicle density.

Exogenous Gonadotropin Rescue of Follicle and Transcriptional Dynamics in Affected Ovaries
In light of the reduced follicle densities and transcriptomic shifts observed in AP and WO-E 2 alligators, we investigated whether these alterations could be attenuated by FSH treatment. To address this, we administered FSH to WO controls, AP, and WO-E 2 animals and examined its ability to rescue follicular (AP and WO-E 2 ) and transcriptomic (AP only) phenotypes ( Figure 4A). Contrary to expectations, reduced follicle densities persisted in AP animals ( Figure 4B-D), and, further, FSH administration did not increase average follicle diameter ( Figure S3) nor density of Stage III follicles in any treatment group compared with non-FSH-treated animals ( Figure S5). However, FSH administration elicited broad transcriptional responses in the ovary ( Figure 4E) and resulted in a partial convergence of AP and WO transcriptomes. Of the 4,169 genes that were initially observed to be highly differentially expressed [log 2 FC >j1j] between AP and WO, FSH administration in AP animals was sufficient to rescue (i.e., reduce the number of DEGs between populations) the expression of 849 genes (20.4%; Figure 4F; white bars, AP+, WO-; numerical figure summary in Table S8). Further, when administered to both AP and WO animals, the degree of FSH-mediated convergence increased, rescuing 2,159 (51.8%) of highly affected genes (white bars, AP+, WO+). Interestingly, the convergence elicited by FSH appeared to be driven predominantly by transcriptional suppression, rather than induction, of affected genes. FSH administration in AP animals rescued fewer than 3% of genes with significantly lower expression in AP animals compared with WO animals (40 of 1,766; blue bar, AP+, WO-), whereas 33.1% of genes with significantly lower expression in AP animals compared with WO were rescued with FSH treatment of both WO and AP animals (blue bar, AP+, WO+). In both instances (AP-FSH compared with WO, and AP-FSH compared with WO-FSH), the proportion of genes with significantly lower expression in AP animals rescued by FSH was smaller than the proportion of genes rescued that were initially elevated at AP. For example, FSH administered to AP animals alone was capable of rescuing 33.7% (809 of 2,403) of highly elevated genes in AP animals compared with WO (black bar, AP+, WO-), and FSH treatment of both AP and WO animals together rescued over 65% of these genes (1,575 of 2,403) (black bar, AP+, WO+).
Consistent with its ability to drive transcriptomic convergence, FSH treatment both in AP animals alone and simultaneous treatment of AP and WO animals together was sufficient to rescue a subset of functional pathways disrupted in AP animals. Of the 62 GO terms initially enriched in genes with significantly higher expression in AP animals, treatment of AP and WO animals with FSH rescued 17 of these affected pathways ( Figure 5A, a and b overlap; Excel Tables S3 and S4), and 2 additional pathways were also rescued by treatment of AP animals alone ( Figure  5A, a, b, and c overlap). Similarly, FSH was capable of rescuing genes and associated pathways with lower expression in AP alligators ( Figure 5B). Of 16 terms enriched in the basal AP-to-WO comparisons, four pathways were rescued by simultaneous treatment of both populations with FSH ( Figure 5B, a and b overlap; Excel Table S5). Interestingly, treatment of AP animals with FSH, but not treatment of both AP and WO animals, was capable A C E D B Figure 2. Assessment of transcriptional profiles in 17b-estradiol (E 2 )-treated juvenile alligators. (A) To explore possible mechanisms underlying populationlevel transcriptomic divergence, Woodruff (WO) embryos were exposed to E 2 prior to the onset of gonadal steroidogenesis (red dot), raised to the juvenile stage under laboratory conditions, and RNAseq against ovarian tissue was conducted (yellow star). (B) Estrogen receptors ESR1 and ESR2 are expressed at the time of dosing, Stage 19, and throughout later stages (20 and 26) of alligator development (n = 4 per stage; bars represent mean ± 1 SD). (C) Volcano plot of differential expression analysis between E 2 -exposed animals and controls; green and blue points are significantly differentially expressed genes between populations (FDR p < 0:05). Green points (within dashed boxes) are significant and are at least 2-fold differentially expressed [log 2 ðFCÞ > j1j] between exposed and control animals. (D) Overlap between differentially expressed genes (DEGs) elevated or suppressed in AP and E 2 -exposed individuals relative to WO controls (WO-VEH). Asterisks denote significance (hypergeometric probability: *** p < 0:0000). (E) Enrichment of ERE-proximate genes in DEGs [log 2 ðFCÞ > j1j] that are suppressed or elevated in AP and in E 2 -exposed animals relative to WO controls (WO-VEH). ERE distance denotes the maximum allowable distance separating a predicted ERE and a gene. Asterisks denote significantly different proportions relative to all detectable ovarian genes (Fisher's exact test: **** p ≤ 0:0001; *** p ≤ 0:001; ** p ≤ 0:01; * p ≤ 0:05). Note: AP, Lake Apopka; CPM, count per million; ERE, estrogen response element; FC, fold change; ESR1, estrogen receptor 1; ESR2, estrogen receptor 2; FDR, false discovery rate; FPT, female-promoting temperature; RNAseq, RNA sequencing; SD, standard deviation; VEH, vehicle. of further rescuing cell cycle and DNA replication pathways (b and c overlap; Excel Table S6), suggesting that FSH may elicit different responses across populations in transcription of genes related to proliferation or the cell cycle. Finally, genes rescued by simultaneous treatment of AP and WO animals with FSH were uniquely enriched for pathways associated with the regulation of lipid biosynthesis and ovarian steroidogenesis (Excel Table S5).

RNAseq qPCR Confirmation
qPCR was employed to assess expression of a subset of genes described herein to independently confirm patterns observed in RNAseq analyses. The genes assessed via qPCR were CYP19A1 ( Figure S6A Figure S6L). RNAseq counts and relative expression values were significantly correlated for 8 of 12 genes (Table S9). Significant differences in expression between FSH-treated and VEH juveniles were detected via two-way ANOVA in RNAseq data for CYP19A1, FST, FSHR, ESR2, AR, AHR1A, AHR2, and GR. These findings were confirmed statistically via qPCR for 4 of 8 genes, CYP19A1, FST, FSHR, and ESR2, and similar effects of FSH were observed for 3 remaining genes, AR, AHR2, and GR expression (although they did not reach statistical significance). Significant expression changes in FSH-treated juveniles were observed in qPCR data for expression of AHR1B, AMH, and PGR that was not observed for RNAseq counts. Similar trends were observed in RNAseq data for these genes that failed to reach significance, however. Significant differences across populations (basal differences, WO-VEH vs. AP-VEH) were detected via unpaired t-test in RNAseq counts for CYP19A1, FST, AHR2, ESR1, AHR1B, FSHR, AR, PGR, and GR; these patterns were confirmed for AHR1B and AHR2 via qPCR. Last, we observed general agreement between RNAseq read counts and qPCR for differences between WO-FSH and AP-FSH-treated juveniles in expression CYP19A1, FST, AHR2, ESR1, ESR2, AR, FSHR, and PGR. These differences were significant in qPCR data for AHR2 and ESR2 and approached significance for FST, ESR1, and FSHR.

Discussion
By employing a nontargeted, transcriptomics-based approach, we a) revealed broad and persistent transcriptional divergence in ovaries from alligators originating from an environment conveying developmental contaminant exposure; b) uncovered functional pathways and signatures of epigenetic modifiers within these patterns of divergence; and c) provided empirical support for a mechanism underlying reproductive pathologies in alligators from AP, wherein exposure to EDCs at a bipotential stage of gonadal development might precociously activate estrogen signaling. Few studies to date have employed transcriptomics to investigate organismal responses to contaminants in environmentally relevant systems (i.e., those composed of genetically diverse individuals in natural populations and exposed to environmental contaminants) (but see Defo et al. 2018;von Hippel et al. 2018), and fewer still have empirically addressed the effects of developmental exposures that persist into later life stages (but see Lea et al. 2016). Nonetheless, the proclivity of environmental endocrine disruptors to shape future ovarian function is well supported . (E) Stage III follicle densities in WO controls, AP, and E 2 -treated WO animals (WO + E2). Letters above box plots denote statistical significance as assessed via one-way ANOVA and Tukey's multiple comparisons test (mean ± 1 SD). Frequency distribution of Pearson correlation coefficients between SIII follicle densities per total cortical area and expression values are depicted for (F) both highly suppressed and highly elevated genes in AP and (G) E 2 -exposed WO animals. Note: ANOVA, analysis of variance; SD, standard deviation.
(for reviews, see Crain et al. 2008;Mahoney and Padmanabhan 2010;Mossa et al. 2017;Padmanabhan et al. 2016) and hinges upon the ability of EDCs to disrupt pathways controlling differentiation and proliferation in the embryonic gonad.
A single dose of E 2 was sufficient to recapitulate AP phenotypes in reference animals, suggesting that the timing of estrogen signaling during development is critical for normal ovarian development. As evidenced by the ontogeny of aromatase expression, gonadal estrogen biosynthesis does not commence until the end of sex determination, approximately four stages after the timing of E 2 exposure used herein (Gabriel et al. 2001;Parrott et al. 2014;Smith et al. 1995). Thus, it is likely that the early gonad develops in the absence of high estrogen levels but remains capable of responding to estrogenic cues. The observed reduction in follicle densities in ovaries from AP alligators suggests that activation of estrogen signaling might be functionally linked to germ cell and/or somatic support cell behavior in the differentiating alligator gonad, which is consistent with evidence that exposure to estrogenic compounds during development can disrupt follicle development (Iguchi 1985;Iguchi and Takasugi 1986;Kim et al. 2009;Stoker et al. 2008). Prior investigations of reproductive abnormalities at AP have uncovered increased prevalence of multi-oocytic follicles (Guillette et al. 1994), an abnormal follicular state linked to precocious estrogenmediated induction of the transforming growth factor-b-signaling components follistatin and inhibin-a. These factors regulate granulosa cell proliferation, germ cell nest breakdown, and follicle Stage III follicle densities in WO controls, AP, and 17b-estradiol (E 2 )-treated WO animals (WO + E2) with exogenous FSH treatment. Letters above box plots denote statistical significance as assessed via one-way ANOVA and Tukey's multiple comparisons test (mean ± 1 SD). (E) Principal components analysis of FSH-treated and nontreated ovarian transcriptomes from WO and AP. (F) Proportion of genes that were initially highly suppressed (blue; 1,766 genes), highly elevated (black; 2,403 genes), or all affected genes (white; 4,169 genes total) in AP animals relative to WO controls that were rescued by exogenous FSH administration; ± indicates which population received FSH to achieve the corresponding proportion of rescued genes. Rescue is defined as a reduction in the number of differentially expressed genes initially detected between groups (i.e., WO-VEH vs. AP-VEH) that was observed when one or more groups were treated with FSH. Asterisks denote significantly different proportions of rescued genes (Fisher's exact test: **** p ≤ 0:0001). Note: ANOVA, analysis of variance; DEGs, differentially expressed genes; FPT, female-producing temperature; PC, principal components; RNAseq, RNA sequencing; SD, standard deviation; VEH, vehicle. assembly during the peri-hatching period (Guillette and Moore 2006). However, this is the first report of reduced Stage III follicles at AP and, taken together with the developmental timing of E 2 treatment, suggests that estrogen might influence follicle development at an earlier stage than previously hypothesized by limiting germ cell proliferation or compromising growth or differentiation of oocytes. Limited evidence exists describing primordial germ cell (PGC) behavior in the alligator, but work by Smith and Joss (1993) has indicated the presence of PGCs in the cortex of the undifferentiated gonad as early as Stage 20, one stage after experimental estrogen dosing in the present study. Furthermore, ERs are actively expressed in the early gonad in alligators (described herein and in Rice et al. 2017), birds (Andrews et al. 1997;Ge et al. 2012;Guioli and Lovell-Badge 2007), turtle (Bergeron et al. 1998;Crews 2007, 2009), eutherian mammals (La Sala et al. 2010), and marsupials (Calatayud et al. 2010), and mammalian data supports receptor expression by PGCs themselves (Calatayud et al. 2010;La Sala et al. 2010). This would collectively suggest that PGCs are present in the gonad during precocious E 2 treatment and are capable of responding to receptor activation.
Although it is true that the study design employed herein is not sufficient to disentangle transcriptional shifts resulting from changes in cellular composition vs. those resulting from cell-autonomous, epigenomic effects of precocious estrogen, consistent enrichment of polycomb components as potential regulators of genes suppressed both in AP and in WO-E 2 animals suggests that epigenetic mechanisms contribute to EDC-induced reproductive pathologies in the alligator. Precocious estrogen signaling could represent a mistimed regulatory cue that disrupts the ontogeny of epigenetic patterning during gonadal development. Consistent with this hypothesis, estrogen signaling can induce expression of EZH2, the H3K27 histone methyltransferase component of the PRC2, via genomic ER signaling while also suppressing EZH2 activity through nongenomic signaling (Bhan et al. 2014;Bredfeldt et al. 2010;Greathouse et al. 2012). Critically, this regulatory function is not limited to endogenous estrogens and can be similarly engaged by nonsteroidal estrogens, such as diethylstilbestrol, bisphenol A, and genistein (Bhan et al. 2014;Bredfeldt et al. 2010;Greathouse et al. 2012). Enrichment of PRC components raises intriguing questions regarding the endogenous function and ontogeny of epigenetic regulators during development, particularly in the context of endocrine regulation. A recent study in a species of turtle with temperature-dependent sex determination revealed that the most proximal factor responsive to temperature in the bipotential gonad is KDM6B, the regulatory counterpart to EZH2 that removes H3K27 methylation to activate gene expression (Ge et al. 2018). In contrast to EZH2 however, expression of KDM6B is suppressed by estrogen, which would suggest that estrogen signaling might represent a crucial hinge in the balance of suppressive vs. permissive epigenetic modifications. Although additional work is necessary to disentangle the developmental role(s) for EREs, estrogen signaling, and the epigenome A B Figure 5. Impacts of exogenous follicle-stimulating hormone (FSH) administration on divergent functional pathways. The ability of exogenous FSH treatment in Apopka (AP) animals alone (circles labeled c; AP + FSH, WO-VEH) or simultaneous treatment of AP and Woodruff (WO) animals (circles labeled a; AP + FSH, WO + FSH) to rescue functional pathways (GO:BP, KEGG) enriched in the divergent AP transcriptome was assessed for both (A) highly elevated and (B) suppressed gene lists in initial AP-VEH vs. WO-VEH contrasts. Terms initially enriched in basal contrasts (circles labeled b) absent FSH treatment were compared with enriched pathways in differentially expressed genes (DEGs) rescued by FSH, with overlap between circles representing pathways that are rescued with FSH. Terms not shared between contrasts represent uniquely enriched pathways, or pathways containing genes that are not rescued by FSH. Note: BP, biological process; GO, Gene Ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes; PKA, protein kinase A; VEH, vehicle. from contemporary ER-ERE mediated effects (e.g., altered circulating estrogen levels in juveniles), these findings nonetheless highlight multiple possible mechanisms for precocious estrogen signaling to elicit persistent functional shifts in the ovary.
The pituitary-derived gonadotropin hormone FSH is a critical regulatory factor controlling follicle development and transcription within the ovary. Specifically, FSH is frequently used to treat female infertility in clinical contexts (Dale et al. 1989;Gemzell 1966;Homburg and Howles 1999;White et al. 2018) because it stimulates growth and maturation of ovarian follicles as well as steroid hormone production by inducing expression of steroidogenic enzymes. Treatment with an exogenous gonadotropin herein partially rescued ovarian transcriptional programs in animals originating from AP. Considering the close association of DEGs with follicle densities and the inability of FSH to increase the number of Stage III follicles, these observations demonstrate a surprising uncoupling of transcriptional activity and follicle dynamics in the juvenile alligator ovary. However, this uncoupling is consistent with a prior study showing that in sexually immature individuals, FSH is capable of stimulating transcription of canonical targets (i.e., CYP19A1) but cannot promote follicular growth beyond early (preantral) stages (François et al. 2017). Because FSH treatment resulted in incomplete transcriptomic convergence in AP animals, broad transcriptional shifts associated with population are not likely driven by outright hypogonadotropism or ovarian insensitivity to FSH. Collectively, these findings suggest that although FSH is capable of reconciling transcriptional patterns across populations, it is not sufficient to rescue reduced follicle densities in AP animals. Thus, phenotypes observed in AP alligators and, by extension, individuals exposed to precocious estrogens, are not solely due to reduced activation of FSH signaling in the ovary. Together, this illustrates additional avenues of research that are needed to fully characterize the etiology and mechanisms of contaminant-induced pathologies at AP.
We provide an experimental demonstration of a potential mechanism underlying reproductive abnormalities observed in alligators from AP, a site historically characterized with high levels of endocrine-disrupting OCPs. Given that OCPs were not quantified in the present study, additional work is needed to establish present-day concentrations in alligator eggs at AP and to show that they are acting directly through the ER (e.g., partial inhibition of ER activation to ameliorate effects). Furthermore, the detection of a minority of differentially expressed genes between AP and WO reference animals that were not recapitulated with precocious estrogen treatment reveals that nonestrogenic actions could also contribute to reproductive abnormalities at AP. Further, the failure of estrogen to perfectly recapitulate AP patterns of expression could highlight contaminant-mediated responses to selection. Nonetheless, these findings collectively support a largely estrogenic mode of action for contaminants to act as drivers of persistent reproductive perturbations in a natural population. Given the conservation of endocrine function across vertebrate taxa, these results should broadly inform our understanding of the developmental origins of reproductive function and the threats posed to environmental and organismal health by EDCs.