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

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 (Galligan et al. 2019; 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 DDT-metabolite–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, 1999; Heinz et al. 1991; Rauschenberger et al. 2007, 2009; Woodward 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. 2005a, 2008; Moore et al. 2010b, 2010c; Pickford 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\mathrm{\mu }\mathrm{g}/\mathrm{g}$ (egg weight) $\text{estradiol-}17\mathrm{\beta }$ (${\mathrm{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 $\mathrm{AP}\text{-}{\mathrm{E}}_2\text{-treated embryos}$ were used). Stage 19 of alligator development precedes the onset of aromatase expression (Parrott et al. 2014), 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 (Hale et al. 2019). Briefly, animals weighing $<92\mathrm{g}$ were fed daily, animals weighing $92–132\mathrm{g}$ were fed three times per week, and animals weighing $>132\mathrm{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\mathrm{\mu }\mathrm{U}/\mathrm{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\mathrm{mg}/\mathrm{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°\mathrm{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 ($\sim 10\mathrm{mg}$ per sample) were manually homogenized in an AGPC lysis solution (water-saturated acidic phenol, $2\mathrm{M}$ guanidinium, thiocyanate, $95\text{\hspace{0.17em}mM}$ sodium acetate, $12\text{\hspace{0.17em}mM}$ sodium citrate, 0.24% $N$-lauroyl sarcosine, $14.4\mathrm{M}$ beta-mercaptoethanol) with a Retsch ball mill and stainless steel beads at 30 Hz for 2 min. Following phase separation with $0.2\mathrm{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 diethylpyrocarbonate-treated 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.

Alligator RNAseq library preparation and sequencing was conducted by the Georgia Genomics and Bioinformatics Core at the University of Georgia. Briefly, $1–2\mathrm{\mu }\mathrm{g}$ 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) $\ge 7.80$, $\text{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 $\mathrm{WO}\text{-}{\mathrm{E}}_2\text{-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 ($\text{parameters mode}=\text{Union}$, $\text{fragments}=\text{TRUE}$, $\text{singleEnd}=\text{FALSE}$, $\text{ignore}.\text{strand}=\text{FALSE}$). This generated a count matrix with 24,848 unique genes in both sequencing runs.

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\mu \mathrm{m}$ 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\times \text{and\hspace{0.17em}}20\times$ 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.

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 ${\mathrm{E}}_2\text{-treated animals}$), low-expression genes [count per million (CPM) $<1$, $\text{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 $\mathrm{WO}\text{-}{\mathrm{E}}_2\text{\hspace{0.17em}animals}$, library sizes were approximately doubled, thus filtering was relaxed ($\text{CPM}<0.5$, $\text{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 ($\text{BCV}=0.30$) to 0.059 ($\text{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, $\text{robust}=\text{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 $\mathrm{WO}\text{-}{\mathrm{E}}_2\text{-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 ${\mathrm{E}}_2\text{-exposed}$ animals was determined using the hypergeometric distribution ($\mathrm{\alpha }=0.05$). Because total features passing low-expression filters differed slightly between experiments ($n=\mathrm{18,435}$ vs. $n=\mathrm{19,230}$), only shared features passing filtering in both experiments were used in the recapitulation analyses ($n=\mathrm{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 ($\mathrm{\alpha }=0.05$). Residual normality was confirmed via the Shapiro-Wilk test of normality ($\mathrm{\alpha }=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 ($\mathrm{\alpha }=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).

Patterns in RNAseq read count were independently compared via quantitative real-time polymerase chain reaction (qPCR) using matched samples for 12 genes described as part of a preceding study (Hale et al. 2019). Genes assessed via qPCR were aromatase (CYP19A1, LOC102566432), follistatin (FST), ER$\mathrm{\alpha }$ (ESR1), ER$\mathrm{\beta }$ (ESR2), FSH receptor (FSHR), androgen receptor (AR), aryl hydrocarbon receptor 1A (AHR1A), aryl hydrocarbon receptor 1B (AHR1B, LOC102568094), aryl hydrocarbon receptor 2 (AHR2, LOC102576203), progesterone receptor (PGR), glucocorticoid receptor (GR or NR3C1), and anti-Mullerian hormone (AMH). Expression levels were determined using absolute quantification (C1000 thermal cycler CFX96 real-time detection system (BioRad)). Briefly, PCR reactions were conducted using a homebrew SYBR® reaction [$0.2\mathrm{\mu }\mathrm{M}$ forward/reverse primer mix, $50\text{\hspace{0.17em}mM}$ potassium chloride, $20\text{\hspace{0.17em}mM}$ Tris-hydrochloride, 0.5% glycerol, 0.5% Tween®-20, 4% dimethyl sulfoxide, $3\text{mM}$ magnesium chloride, $20\mathrm{\mu }\mathrm{M}$ deoxynucleotide mix, $0.5\times$ SYBR® Green (Life Technologies), and $0.1\mathrm{U}/\mathrm{\mu }\mathrm{L}$ AmpliTaq Gold (Applied Biosystems)] and $2\text{-}\mathrm{\mu }\mathrm{L}$ complementary DNA template. Reactions were prepared in a total volume of $50\mathrm{\mu }\mathrm{L}$, split, and run in $15\text{-}\mathrm{\mu }\mathrm{L}$ triplicates. Expression levels were determined by interpolation of a plate-matched standard curve of plasmid controls of known concentration for each gene and subsequently normalized to the geometric mean of ribosomal protein L8 (RPL8) and eukaryotic elongation factor 1 (EEF1) expression levels. Any triplicate exhibiting coefficient of variation levels $>0.4$ prior to normalization was excluded from analyses. All primers were intron spanning and sequences, along with annealing temperatures and amplicon lengths, are reported in Table S2.

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) ($\mathrm{\alpha }=0.05$). Significance of effects were then compared across quantification approaches. Residual normality and homoscedasticity were confirmed via the Shapiro-Wilk test ($\mathrm{\alpha }=0.05$) and Spearman’s test ($\mathrm{\alpha }=0.05$), respectively, for both approaches. Values were transformed as necessary to meet these model assumptions.

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, 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) ($\mathrm{e}\text{-value cutoff}=1\times {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 ${\mathrm{E}}_2$ experiment. Any LOC loci that was not identified via this approach was excluded from downstream analyses.

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\times {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\text{\hspace{0.17em}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 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\text{\hspace{0.17em}fold change\hspace{0.17em}}\left(\text{FC}\right)>|1|$] was compared with total enrichment in all annotated ovarian genes passing filtering (background) using Fisher’s exact tests ($\mathrm{\alpha }=0.05$) [Table S3 (AP-responsive) and Table S5 (${\mathrm{E}}_2\text{-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.

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 ($\mathrm{\alpha }=0.05$).

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\mathrm{d}$), and 26 (Stage $19+30\mathrm{d}$) incubated at a female-promoting 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.

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 ${\mathrm{E}}_2$ at Stage 19 (Figure 2A), a stage that precedes the detection of aromatase expression and sexually dimorphic gonadal steroidogenesis (Kohno et al. 2014; 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-${\mathrm{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\text{-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 ${\mathrm{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 ${\mathrm{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\text{FC}>|1|$] in varying degrees of proximity to predicted EREs in the alligator genome. We found that genes with significantly greater expression in AP or ${\mathrm{E}}_2\text{-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 $\mathrm{WO}\text{-}{\mathrm{E}}_2\text{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).

We next sought to investigate follicular profiles in AP and ${\mathrm{E}}_2\text{-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_{\mathrm{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 ${\mathrm{E}}_2$ was sufficient to phenocopy the lowered follicle density observed in animals from AP given that both AP and ${\mathrm{E}}_2\text{-treated}$ WO animals ($\mathrm{WO}\text{-}{\mathrm{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. $\mathrm{WO}\text{-}{\mathrm{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\left(\text{FC}\right)<–1$] in AP and $\mathrm{WO}\text{-}{\mathrm{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\left(\text{FC}\right)>1$] were positively associated with follicle density.

In light of the reduced follicle densities and transcriptomic shifts observed in AP and $\mathrm{WO}\text{-}{\mathrm{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 $\mathrm{WO}\text{-}{\mathrm{E}}_2$ animals and examined its ability to rescue follicular (AP and $\mathrm{WO}\text{-}{\mathrm{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\text{FC}>|1|$] 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, $\text{AP}+$, $\text{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, $\text{AP}+$, $\text{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, $\text{AP}+$, $\text{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, $\text{AP}+$, $\text{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, $\text{AP}+$, $\text{WO}–$), and FSH treatment of both AP and WO animals together rescued over 65% of these genes (1,575 of 2,403) (black bar, $\text{AP}+$, $\text{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, $\mathrm{\alpha }$ and $\mathrm{\beta }$ overlap; Excel Tables S3 and S4), and 2 additional pathways were also rescued by treatment of AP animals alone (Figure 5A, $\mathrm{\alpha }$, $\mathrm{\beta }$, and $\mathrm{\gamma }$ 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, $\mathrm{\alpha }$ and $\mathrm{\beta }$ overlap; Excel Table S5). Interestingly, treatment of AP animals with FSH, but not treatment of both AP and WO animals, was capable of further rescuing cell cycle and DNA replication pathways ($\mathrm{\beta }$ and $\mathrm{\gamma }$ 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).

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), FST (Figure S6B), ER$\mathrm{\alpha }$ (Figure S6C; ESR1), ER$\mathrm{\beta }$ (Figure S6D; ESR2), FSHR (Figure S6E), AR (Figure S6F), AHR1A (Figure S6G), AHR1B (Figure S6H), AHR2 (Figure S6I), PGR (Figure S6J), GR (Figure S6K), and AMH (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.