Genome‐wide DNA methylation patterns harbour signatures of hatchling sex and past incubation temperature in a species with environmental sex determination

Abstract Conservation of thermally sensitive species depends on monitoring organismal and population‐level responses to environmental change in real time. Epigenetic processes are increasingly recognized as key integrators of environmental conditions into developmentally plastic responses, and attendant epigenomic data sets hold potential for revealing cryptic phenotypes relevant to conservation efforts. Here, we demonstrate the utility of genome‐wide DNA methylation (DNAm) patterns in the face of climate change for a group of especially vulnerable species, those with temperature‐dependent sex determination (TSD). Due to their reliance on thermal cues during development to determine sexual fate, contemporary shifts in temperature are predicted to skew offspring sex ratios and ultimately destabilize sensitive populations. Using reduced‐representation bisulphite sequencing, we profiled the DNA methylome in blood cells of hatchling American alligators (Alligator mississippiensis), a TSD species lacking reliable markers of sexual dimorphism in early life stages. We identified 120 sex‐associated differentially methylated cytosines (DMCs; FDR < 0.1) in hatchlings incubated under a range of temperatures, as well as 707 unique temperature‐associated DMCs. We further developed DNAm‐based models capable of predicting hatchling sex with 100% accuracy (in 20 training samples and four test samples) and past incubation temperature with a mean absolute error of 1.2°C (in four test samples) based on the methylation status of 20 and 24 loci, respectively. Though largely independent of epigenomic patterning occurring in the embryonic gonad during TSD, DNAm patterns in blood cells may serve as nonlethal markers of hatchling sex and past incubation conditions in conservation applications. These findings also raise intriguing questions regarding tissue‐specific epigenomic patterning in the context of developmental plasticity.


| INTRODUC TI ON
Climate change threatens species persistence in wild populations via its influence on key demographic parameters. In ectotherms, shifting temperature regimes impact juvenile recruitment (Laloë et al., 2017;Noble et al., 2018;Pruett & Warner, 2021;Schwanz et al., 2010), adult fecundity (Crews et al., 1998;Wang et al., 2016;Warner & Shine, 2008) and, in some cases, population sex ratios (Bock et al., 2020;Janzen, 1994;Jensen et al., 2018;Reneker & Kamel, 2016). Fortunately, efforts to predict organismal and population-level responses to changing environmental conditions are rapidly advancing with improvements in the precision and geographical resolution of global climate models and their integration with empirical data from physiological, ecological and evolutionary studies (Carlo et al., 2018;Laloë et al., 2014;Mitchell et al., 2008;Riddell et al., 2018;Seebacher et al., 2015). Yet, one major challenge moving forward will be ground truthing these predictions to adapt conservation strategies accordingly. Genomic and epigenomic data sets represent powerful resources for conservationists in these efforts (Formenti et al., 2022;Layton & Bradbury, 2022;Waldvogel et al., 2020). The use of genomic tools in conservation is already widespread, allowing for assessments of population structure and connectivity (Wright et al., 2020), inbreeding depression (Dussex et al., 2021;Kardos et al., 2016), contemporary evolutionary responses to environmental change (Bi et al., 2019;Catullo et al., 2019) and adaptive capacities to cope with future change (Bay et al., 2018;Flanagan et al., 2018;Harrisson et al., 2014;Layton et al., 2021). The potential of epigenomic tools, on the other hand, is only just beginning to be realized. DNA methylation (DNAm), the covalent addition of a methyl group to a cytosine base, is one of the best-studied epigenetic marks, and variation in DNAm within and across individuals harbours key biological and ecological insights (Angers et al., 2010;Schübeler, 2015;Ziller et al., 2013). In the face of a rapidly shifting global climate, variation in genomic DNAm patterns has the potential to reveal cryptic organismal traits relevant to conservation and may provide insight into the environmental conditions from which they arose.
The potential of DNAm patterns to serve as biomarkers for environmentally induced organismal phenotypes is particularly relevant for thermally sensitive species facing climate-induced population collapse. Among these vulnerable taxa are species with temperature-dependent sex determination (TSD). Due to their reliance on temperature cues during incubation to determine offspring sex, species with TSD, including many turtles, all crocodilians and some fish, are predicted to experience highly skewed population sex ratios as a result of climate change (Bock et al., 2020;Janzen, 1994;Laloë et al., 2014;Mitchell et al., 2008Mitchell et al., , 2010Ospina-Alvarez & Piferrer, 2008;Reneker & Kamel, 2016). Such skews in sex ratios are already being detected in wild populations of TSD species (Jensen et al., 2018) and may exert lasting adverse impacts on population dynamics given the long generation times of many TSD species (Böhm et al., 2016). Despite these observations, efforts to actively monitor cohort sex ratios in the field remain limited. This is in part due to the lack of reliable, nonlethal methods of sexing individuals at early life stages. Since hatchlings do not tend to exhibit sexually dimorphic morphology, specialized training is required to determine hatchling sex (Allsteadt & Lang, 1995), and genetic markers of sex (i.e., sexlinked genes or alleles) are nonexistent in TSD species (or may not be indicative of phenotypic sex in the case of temperature-induced sex reversal; Holleley et al., 2015). DNAm patterns in blood represent a novel avenue by which those tasked with managing populations of TSD species under climate change might gain insight into ongoing sex ratio shifts as well as changes in field incubation conditions.
The DNA methylome appears to integrate temperature cues into sexually dimorphic phenotypic trajectories during TSD, though the details of this process have yet to be resolved at a genome-wide scale. During development, gonadal tissues acquire sex-associated DNAm patterns at the promoters of genes with conserved roles in vertebrate sex determination in response to incubation temperature (Matsumoto et al., 2016;Navarro-Martín et al., 2011;. It remains unknown whether somatic tissues that lack direct roles in sex determination also acquire similar sex-or temperature-associated patterns. However, a few lines of evidence from species with genotypic sex determination support this possibility. Several studies have identified autosomal sex-associated variation in DNAm patterns across species differing in sex-determination system in somatic tissues (Gatev et al., 2021), including blood (Caizergues et al., 2021;Janowitz Koch et al., 2016), muscle (Wan et al., 2016) and liver (Hu et al., 2019;Zhuang et al., 2020). Some sex-biased autosomal DNAm patterns have even been shown to be established during development and stable throughout the life course (Gatev et al., 2021). Though such differences in DNAm patterns tend to be more subtle in autosomal regions of the genome compared to sex chromosomes, these findings suggest the DNA methylomes of somatic tissues may harbour sex-associated patterns in TSD species. Further, environmental influences on DNAm patterns established during development can persist in the postnatal period even after the environmental cue has passed (Gavery et al., 2019;von Holdt et al., 2021). Developmental thermal cues can exert lasting influences on the DNA methylome (Anastasiadi et al., 2021;Metzger & Schulte, 2017), raising the possibility that DNAm patterns in hatchlings might provide information about past incubation conditions.
In this study, we investigate sex-and temperature-associated variation in the DNA methylomes of hatchling blood cells in the American alligator, a species with TSD. DNAm dynamics in blood cells are examined in the context of other sex-linked biological correlates such as circulating levels of sex steroids and gonadal DNAm patterns established during TSD. We further use these DNAm patterns to develop models that predict hatchling sex and past incubation temperature based on the methylation status of a subset of loci. Using these predictive models, we assess the potential utility of DNAm-based markers in the conservation of species with TSD under rapid environmental change. Due to the unique temperatureby-sex ratio reaction norm of the alligator, this species serves as a particularly insightful model to identify sex-and temperatureassociated DNAm patterns. As in all other crocodilians, females are produced at both low and high incubation temperatures, while males are produced at intermediate incubation temperatures (Lang & Andrews, 1994). From the perspective of experimental design, this pattern of TSD allows us to disentangle DNAm variation linked to sex from that linked to temperature through the inclusion of females from both the low and high ranges of the reaction norm. The shape of the crocodilian temperature-by-sex ratio reaction norm also yields greater uncertainty in predictions of sex ratio shifts in response to climate change (Bock et al., 2020). Thus, monitoring sex ratios in the field is especially pertinent for these species and our results hold unique relevance for future conservation.

| Incubation experiment and sample collection
In June 2013, seven clutches of Alligator mississipiensis eggs were collected from Lake Apopka, Florida, and transported to Hollings Marine Laboratory in Charleston, South Carolina (McCoy et al., 2016). Detailed methods regarding the incubation experiment can be found in the associated report, which sought to decipher temperature-dependent inter-and intrasexual variation in gonadal gene expression . Briefly, eggs were kept in damp sphagnum moss within commercial incubators (Thermo Scientific, Forma Environmental Chambers, model #3920). One embryo from each clutch was staged according to Ferguson developmental staging criteria (Ferguson, 1985). Viable eggs were maintained at 30°C until stage 19, the opening of the canonical thermosensitive period (Lang & Andrews, 1994), after which eggs from each clutch were distributed across four incubation temperatures (30, 32, 33.5, 34.5°C), where they were maintained until hatch. All animal protocols were approved by the Medical University of South Carolina's Institutional Animal Care and Use Committee (AUP # 3036).
Following hatch, alligator neonates were measured to assess morphometric traits including mass, snout-vent length and total length. At 7 days post-hatch, body morphometric measurements were repeated, and a blood sample was taken from the postcranial sinus. Animals were then euthanized with a lethal dose of pentabarbital (Sigma Aldrich) and necropsied to collect the whole brain and gonad-adrenal-mesonephros (GAM) complex. Fresh tissues were preserved in RNAlater (Life Technologies) and kept at −80°C. Whole blood samples were stored in lithium heparin Vacutainer tubes (BD) on ice prior to centrifugation at 1500 rcf for 20 min. Plasma was isolated and RNAlater (Life Technologies) was added to the residual pelleted blood cells. Plasma was stored at −80°C and blood cells in RNAlater were stored at −20°C. Plasma concentrations of 17βoestradiol (E2) and testosterone (T) were quantified via radioimmunoassays as described in McCoy et al. (2016).

| Temperature-by-sex reaction norm, sex ratio assessment, and sample selection
To place the incubation treatments and resulting hatchling sex ratios of the present study in the broader context of historical incubation experiments, we gathered data from 14 studies in the American alligator employing constant incubation temperatures and quantifying hatchling sex ratios (Table S1). An equation for the temperatureby-sex ratio reaction norm in this species was then fit with these data using logistic regression with a quadratic term for temperature ( Figure 1).
In the present study, hatchling sex was assigned based on gonadal expression of AMH and CYP19A1 in McCoy et al. (2016). These transcripts exhibit robust sex-biased expression patterns in hatchling gonads, and in a previous study, sex assignment based on expression of these genes was confirmed by histology in 98.7% (213/216) of samples (Kohno et al., 2015). Consistent with previous reports of the alligator temperature-by-sex reaction norm (Lang & Andrews, 1994) Figure 1). While incubation temperatures exceeding 34°C are generally regarded as female-promoting, considerable variation is observed in the sex ratios resulting from these high temperatures (Figure 1). While the factors contributing to this variation (e.g., genetic variation, differences in the early incubation environment) are poorly understood, this observation underscored the importance of including hatchlings from high incubation temperatures in the current analyses aimed at developing DNAm-based markers of hatchling sex. To identify sexually dimorphic DNAm patterns in blood cells from hatchlings spanning the entire temperatureby-sex ratio reaction norm, samples from eight females incubated at 30°C, eight males incubated at 33.5°C, four females incubated at 34.5°C and four males incubated at 34.5°C were selected for reduced representation bisulphite sequencing (RRBS; Table 1). Clutch identity was also accounted for in sample selection, ensuring that clutch representation (N = 5 clutches included) was balanced across incubation temperatures and sexes.

| Nucleic acid isolation and quantification
DNA was isolated from blood cells using a modified column approach based on the SV total RNA isolation system (Promega). Frozen samples were thawed on ice. To remove RNAlater from the preserved cells, 500 μl of phosphate-buffered saline (PBS) was added to a 200μl aliquot of each sample and samples were centrifuged at 5000 rcf for 15 min at 4°C (and for an additional 5 min if a cell pellet did not form). Pelleted cells were washed once more with 500 μl of PBS and nucleic acids were isolated according to the protocol reported in Bae et al. (2021); Appendix S1). Resulting DNA was eluted in 60 μl TE buffer (pH 8.0). DNA concentrations and qualities were quantified using a Qubit fluorometer (Thermo Fisher) and NanoDrop spectrometer (Thermo Fisher) respectively. Extracted DNA was stored at −80°C until library preparation.

| Reduced-representation bisulphite sequencing
Reduced representation bisulphite sequencing libraries were constructed with 200 ng input DNA from each sample (n = 24; Table 1

| Quality trimming, read alignment and methylation data processing
Raw reads were trimmed to remove adapter contamination and low-quality bases (Phred score < 20) using trimgalore! (version 0.6.5) with the "--rrbs" option and a minimum overlap with the adapter sequence of 1 bp ("-stringency 1"; https://www.bioin forma tics. babra ham.ac.uk/proje cts/trim_galor e/). After quality trimming, F I G U R E 1 Incubation temperature treatments in the context of the American alligator temperature-by-sex ratio reaction norm. Each point corresponds to the resulting sex ratio from a single incubation treatment from one of 14 experiments. Points corresponding to treatments and resulting sex ratios used in this study are coloured according to their temperature (30°C = LFPT, blue; 33.5°C = MPT, green; 34.5°C = HFPT, yellow). Point size is proportional to the number of hatchlings that were sexed. Solid line indicates the fit of the logistic model of sex ratio with a quadratic term for temperature, and the shaded area indicates the 95% confidence interval.  Resulting BAM files were sorted and indexed with samtools (version 1.10; Li et al., 2009), then read into R (version 4.0.0; R Core Team, 2021) with the "processBismarkAln" function in the methylkit package (version 1.14.2; bioconductor version 2.48.0; Akalin et al., 2012).
Analyses were limited to cytosines in a CpG context, as most cytosine methylation in vertebrate genomes occurs in this context (Feng et al., 2010

| Differential methylation analysis
Differentially methylated cytosines (DMCs) were identified via five main comparisons. To identify "sex-associated" DMCs, all females were compared to all males (FvM) and 34.5°C females were compared to 34.5°C males (F 34.5 vM 34.5 ). If a cytosine was identified as significantly differentially methylated in the FvM contrast or the F 34.5 vM 34.5 contrast, it was considered "sex-associated." Further, a subset of these "sex-associated" sites were considered "universal sex-associated" sites if they were significantly differentially methylated in the FvM contrast (Appendix S1). To identify "temperatureassociated" DMCs, pairwise comparisons of hatchlings from each temperature treatment were conducted (30v33.5, 30v34.5, 33.5v34.5). If a cytosine was identified as differentially methylated in both the sex-and temperature-based differential methylation analyses, the "sex-associated" definition took precedence over the "temperature-associated" definition. By treating all CpG sites identified as differentially methylated in both the sex-and temperaturebased differential methylation analyses as "sex-associated," we ensured that those sites defined as "temperature-associated" only included sites with differential methylation patterns related to temperature, rather than the confounding effect of sex (i.e., a DMC identified in the 30v33.5 contrast could be related to sex or temperature because these two variables are completely confounded, but if the site was also identified as a DMC in the FvM or F 34.5 vM 34.5 comparison, we assume it is related to sex). Differential methylation was modelled in methylkit using logistic regression and overdispersion of variance was corrected for according to McCullagh and Nelder (1989). A cytosine was defined as differentially methylated with respect to the focal comparison if it passed a 10% Benjamini-Hochberg false discovery rate (FDR) threshold. Resulting "sexassociated" DMCs were categorized based on whether percent methylation was greater in females or males, while "temperatureassociated" DMCs were categorized based on whether percent methylation increased or decreased with incubation temperature.
We also assessed temperature-associated DNAm patterns by performing Spearman correlations between incubation temperature and percent methylation for all covered cytosines using the "co-rAndPvalue" function within the wgcna package (version 1.70-3; Langfelder & Horvath, 2008). To correct for multiple comparisons, an FDR approach was applied using the "p.adjust" function in R. were compared to those of loci not identified as temperatureassociated DMCs using a nonparametric Mann-Whitney U test.

| Annotation of differentially methylated loci
Differentially methylated cytosines were annotated according to their gene context (i.e., promoter, exon, intron or intergenic) using the alligator genome annotation (GCF_000281125.3, annotation release 102) within the genomicranges package in R. If a cytosine overlapped with multiple gene contexts, the site was defined according to the following precedent: promoter > exon > intron > intergenic.
The alligator gene annotation table from NCBI was used within genomicranges to determine the closest gene to each cytosine in the data set using the "distancetoNearest" function. Gene ontology (GO) term enrichment analysis was performed for unordered gene lists associated with DMCs identified as temperature-associated, sex-associated and universal sex-associated using g:profiler (version e105_eg52_p16_e84549f; Raudvere et al., 2019). Separate GO analyses were performed for gene lists corresponding to DMCs overlapping promoters and gene lists corresponding to DMCs overlapping promoters or gene bodies. GO terms were considered significantly enriched if they passed a Benjamini-Hochberg FDR threshold of 10% against a custom background composed of genes associated with all analysed cytosines. Any genes lacking annotations (i.e., uncharacterized "LOC" genes) were excluded from analysis.
The regulatory potential and functional role of individual CpG sites is often linked to local CpG density, especially in the context of CpG "islands" in mammals (Weber et al., 2007). Thus, to better understand the relevance of DMCs identified in our analyses, we annotated these loci with respect to CpG density via two complementary approaches. First, we predicted the locations of CpG islands (CGIs) in the alligator genome using the program "cpgplot" (emboss version 6.6.0; parameters: -window 100, -minlen 200, -minoe 0.6, -minpc 50; Rice et al., 2000). Previous work has demonstrated that CpGs with greater overall methylation tend to be more dynamic than those with low overall methylation (Ziller et al., 2013). To examine evidence for this pattern in the present study, DMCs were also characterized according to the overall percent methylation at these loci. Mean percent methylation across all samples was determined for all cytosines analysed, then loci were categorized into one of five bins (mean percent methylation = 0%-20%, 20%-40%, 40%-60%, 60%-80% or 80%-100%). To Appendix S1). Fisher's exact tests followed by FDR correction were conducted for each annotation category: gene context, CGI context, CpG density context and overall per cent methylation context.

| Potential regulation by circulating sex steroids
To investigate the potential contribution of circulating sex steroid hormones towards observed differential methylation patterns, we performed Spearman correlations between percent methylation for all covered cytosines and plasma E2 using the "corAndPvalue" in the wgcna package. We also tested for Spearman correlations between percent methylation and plasma testosterone. Correlations were tested with the imputed percent methylation data set, and correction for multiple comparisons was performed using an FDR approach (FDR < 0.1). The absolute Spearman correlation coefficients with E2 and testosterone for loci identified as sex-associated DMCs were compared to those of all other covered loci using nonparametric Mann-Whitney U tests.
To further assess the potential role of sex steroids in contributing to differential methylation patterns, we examined the proximity of DMCs to putative hormone response elements. The locations of putative oestrogen-response elements (EREs) and androgenresponse elements (AREs) in the alligator genome were predicted using the tool pwmscan, a genome-wide position weight matrix scanner (Ambrosini et al., 2018)

| DNAm-based predictor of sex and incubation temperature
To develop a DNAm-based model to predict hatchling sex, we divided our samples into a training set of 20 individuals and a test set of four individuals. The test set was composed of one randomly selected individual from each temperature-sex combination. Loci exhibiting differential methylation with respect to sex across all individuals in the training set were identified as described previously, but in this case, cytosines were considered differentially methylated if they passed a more stringent FDR threshold of 5%. We then identified the subset of cytosines that were differentially methylated with respect to sex within the training set and across the entire data set. The imputed percent methylation matrix for these loci was then modelled using penalized binomial logistic regression with elastic net regularization (alpha = 0.5) in the R package glmnet (version 4.1-3) to identify sites for which methylation status was most predictive of sex in the training set (Zou & Hastie, 2005). In this modelling approach, the coefficients of loci that are weak predictors of sex are reduced to zero, thereby eliminating them from the model. The strength of the penalty on model coefficients is controlled by the tuning parameter lambda, which is selected to minimize mean-squared error during internal n-fold cross-validation (n-fold = 5). The performance of the model was then assessed by the correct assignment of sex for individuals in the test set. As a point of comparison, the model was also fit using a smaller training set of 12 samples and tested on the other 12 samples (Appendix S1).
To develop a DNAm-based predictive model for incubation temperature, a similar approach was taken to the one described for developing a predictive model of sex. Loci exhibiting differential methylation with respect to temperature were used as input for the predictor selection. In this case, we fit a penalized linear regression with elastic net regularization. Model fit on the training set and performance on the test set was evaluated through the R 2 of the relationship between actual and predicted incubation temperature and the absolute error of the predicted incubation temperatures.

| Comparison of DNAm patterns between blood and embryonic gonad
To place the temperature-associated DNAm patterns of hatchling blood cells in the broader context of global epigenetic patterning during TSD, we utilized an additional data set consisting of RRBS data from alligator embryonic gonads sampled just after the closing of the thermosensitive period in development. With these data, we assessed the extent to which genome-wide DNAm responses to incubation temperature differ between somatic (i.e., blood) and gonadal tissues. Briefly, seven clutches of alligator eggs (N = 315) were collected from wild nests in June 2018 at the Lake Woodruff High Output flowcell (75 bp, single-end). Sequencing generated a total of 251 million reads, with 93% meeting or exceeding the Phred score quality threshold of >30. Raw reads were processed through the same bioinformatic pipeline as described previously. One sample was excluded from further analysis due to aberrant genomic coverage. Temperature-associated DMCs were then identified in the gonad through pairwise comparisons of individuals from each temperature group using logistic regression with correction for overdispersion as described previously (FDR < 0.1).
To assess the extent to which temperature-associated differential methylation patterns observed in the hatchling blood cells overlapped with those in the embryonic gonad, we subset each data set to only include the cytosines that met coverage and filtering requirements in both tissues (N = 181,834). We compared the proportion of shared loci that were differentially methylated with respect to temperature across tissues using a Fisher's exact test. The distribution of temperature-associated DMCs with respect to gene context was also compared across tissues using Fisher's exact tests followed by

| Genome-wide DNAm patterns harbour signatures of hatchling sex and past incubation temperature
A total of 462,236 cytosines were retained in our analyses following filtering (denoted "covered cytosines" or "covered CpGs"). Covered CpGs make up a relatively small proportion (2.2%) of the 21,480,261 CG dinucleotides in the alligator genome. However, as has been previously demonstrated with RRBS techniques, covered CpGs were enriched in regions of potential functional significance relative to all the CpGs in the genome (Appendix S1; Figure S1), including in promoters (p < .0001; Figure S1A) and exons (p < .0001; Figure S1A), as well as CpG islands (p < .0001), CpG shores (p < .0001) and CpG shelves (p < .0001). Of the CpGs analysed, 120 (0.026%) loci dis- Further, universal sex-associated DMCs were enriched in highly methylated loci (80%-100% methylation in all samples; p = .039; Figure S1C). When we relaxed the requirement that loci be differentially methylated across all females and males, we identified 1307 (0.28%) sex-associated DMCs (significant in FvM or F 34.5 vM 34.5 ).
This group included loci that were only differentially methylated between sexes in hatchlings incubated at 34.5°C. Interestingly, only 10 CpGs were identified as differentially methylated across all females and males as well as within 34.5°C individuals, suggesting that a large proportion of sex-associated DMCs (n = 1187; 90.8%) are specific to hatchlings incubated at 34.5°C.
Of the CpGs analysed, 707 (0.15%) loci exhibited differential methylation with respect to incubation temperature and not sex ( Figure 3). Similar proportions of these temperature-associated DMCs showed greater methylation with increasing temperatures (n = 356; 50.4%) and greater methylation with decreasing temperatures (n = 351; 49.6%). The mean difference in percent methylation between incubation temperature groups at these DMCs was 29.9% for loci that exhibited increased methylation with increased temperature and 29.5% for loci that exhibited decreased methylation with increased temperature (Figure 3c). Similar to universal sex-associated DMCs, temperature-associated DMCs tended to be highly methylated and were enriched in loci with an overall mean percent methylation of 80%-100% relative to covered CpGs (p < .0001; Figure S1C). The comparison of 30°C hatchlings to 34.5°C yielded the greatest number of temperature-associated DMCs (n = 273), while the comparison of 30°C hatchlings to 33.5°C hatchlings yielded the least (n = 230). Few loci were identified as differentially methylated with respect to temperature across multiple pairwise temperature comparisons, and none were identified as significant across all three temperature comparisons. When testing for Spearman correlations between temperature and methylation status across all covered CpGs, we did not identify any significant correlations after FDR correction ( Figure 4). Even so, several loci showed robust correlation coefficients with temperature (absolute correlation coefficients of ~0.8; top eight most correlated loci shown in Figure 4c). The absolute correlation coefficients of DMCs identified as temperature-associated through the pairwise comparisons were significantly greater than those of loci not identified as temperatureassociated (p < .0001; Figure 4b).

| DMCs do not show enrichment in promoters, but those that do overlap promoters and gene bodies tend to associate with genes involved in blood-and immune-related processes
If the DMCs we identified are involved in gene regulation, we might expect them to be enriched in genomic regions of functional relevance such as promoters or gene bodies. We did not, however, observe such enrichment. In fact, sex-associated DMCs were F I G U R E 2 Sex-associated differential methylation patterns. (a) Violin plots depicting kernel density of the per cent methylation of all universal sex-associated differentially methylated cytosines (DMCs; n = 120) for each treatment group, with panels separating DMCs that are hypermethylated in females (female-biased) from those that are hypermethylated in males (male-biased). Central points of violins represent overall mean per cent methylation. significantly depleted in promoter regions (p = .00048) relative to covered CpGs. On the other hand, intergenic CpG sites may correspond to cis-regulatory regions (Jones, 2012;Stadler et al., 2011;Stroud et al., 2011;Wiench et al., 2011) and sex-associated DMCs were enriched in intergenic regions (p = .013; Figure S1A). This pattern was not observed when only considering the universal sexassociated DMCs. These DMCs, as well as temperature-associated DMCs, exhibited a similar distribution across different gene contexts to all the covered CpGs ( Figure S1A), with subsets of DMCs overlapping with promoters (9.2% of universal sex-associated DMCs, 6.5% of sex-associated DMCs and 7.6% of temperature associated DMCs) and gene bodies (45.8% of universal sex-associated DMCs, 39.3% of sex-associated DMCs and 41.6% of temperature-associated DMCs).
Genes with sex-associated DMCs within their promoters were enriched for the GO term "DNA binding domain binding" (GO:0050692; Table S3), and genes with temperature-associated DMCs within their promoters were enriched for the GO terms "integral component of synaptic membrane" (GO:0099699), "cell-cell contact zone" (GO: 0044291), "early phagosome" (GO:0032009) and "intrinsic component of synaptic membrane" (GO:0099240; Table S3). When we considered genes with DMCs within their promoters or gene bodies, we identified several more enriched GO terms, with the blood-related GO terms "vasculature development" (GO:0001944) and "blood vessel development" (GO:0001568) enriched for genes with temperature-associated DMCs within their promoters or gene bodies (Table S4). concentrations, though the effect of sex was close to the threshold for significance (α = 0.05; F 1,18.91 = 3.87; p = .064; Figure 5a).

| Sex-associated DMCs collectively exhibit correlations with oestradiol concentrations
We also did not detect significant Spearman correlations between plasma E2 and the methylation status of any one cytosine after the FDR correction, though several loci showed robust correlation coefficients (top eight most correlated loci shown in Figure 5d). Sexassociated DMCs did, however, show significantly greater absolute correlation coefficients with E2 concentrations than loci that were not sex-associated (Figure 5b,c). The same was not observed for T concentrations ( Figure S3). The apparent link to oestrogen signalling was not reflected in the proximity of sex-associated DMCs to putative EREs, nor AREs, as this did not differ significantly from that of covered CpGs (Figures S2 and S4).

| Methylation patterns of separate subsets of cytosines predict hatchling sex and past incubation temperature
Beyond describing sex-and temperature-associated DNAm patterns, we also aimed to determine whether DNAm patterns in blood F I G U R E 3 Temperature-associated differential methylation patterns. (a) Violin plots depicting kernel density of the per cent methylation of all temperature-associated differentially methylated cytosines (DMCs; n = 707) for each treatment group, with panels separating DMCs that exhibit increased methylation with increasing temperature (positive) from those that exhibit decreased methylation with increasing temperature ( (Table S5). This model correctly assigned sex to all 20 of the training set individuals with minimal variability in the model predictions (mean ± SE predicted probability male: females = 0.07 ± 0.01; males = 0.93 ± 0.01; Figure 6a). This model also correctly assigned sex to all four of the test set individuals, though there was more variability in the model predictions (mean ± SE predicted probability male: females = 0.28 ± 0.07; males = 0.71 ± 0.01; Figure 6b). For the predictive model of past incubation temperature, 152 cytosines passed the threshold to be included in variable selection, and 24 of those cytosines were selected to be included in the final model (Table S6).
This model fit the training set well, with an R 2 > .99 between the actual incubation temperatures and fitted incubation temperatures.
The mean absolute error for the training set individuals was 0.05°C ( Figure 6c). However, the model did not perform as well on the test set individuals. The R 2 for the test set was .81 and the mean absolute error of the model predictions was 1.20°C (Figure 6d). Splitting the data into equally sized training and test sets (12 samples in each) did not substantially alter these results (Appendix S1).

| Embryonic gonads following temperaturedependent sex determination show more extensive temperature-associated DNAm patterning compared to hatchling blood cells
Though we observed both sex-and temperature-associated patterns in the DNA methylomes of hatchling blood cells, these patterns were considerably less extensive than the temperatureassociated DNAm patterns we observed in the embryonic gonad following TSD (Figure 7). In our comparison of the blood and gonadal tissues, we used information on 181,834 CpG sites (39.3% of covered CpGs in blood) which overlapped both RRBS data sets (Table S7). Of these shared, covered cytosines, 232 (0.13%) were identified as differentially methylated with respect to incubation temperature in blood, whereas 4339 (2.39%) cytosines were identified as differentially methylated with respect to incubation temperature in the gonad (Figure 7a). The proportion of DMCs identified in the gonad was significantly greater than that identified in blood (p < .0001). Only five loci were identified as temperature-associated DMCs in both tissues and the direction of the methylation difference was not consistent across tissues. The maximum absolute difference in percent methylation for DMCs identified in the gonad was 100%, while the maximum absolute different in percent methylation for DMCs identified in the blood was 66.67% (Figure 7b). Temperature-associated DMCs identified in the blood displayed a similar genomic distribution with respect to gene context as all shared, covered CpGs ( Figure 7c). However, temperature-associated DMCs in the gonad were depleted in promoters and enriched in introns relative to all shared, covered CpGs ( Figure 7c).

| DISCUSS ION
In the present study, we demonstrate that the DNA methylomes of The functional relevance of DNAm for environmentally sensitive transcriptional regulation tends to be highly dependent on the genomic context in which it occurs (Jones, 2012;Schübeler, 2015).
While promoter methylation has traditionally been regarded as a marker of transcriptional repression, further work has revealed more nuanced relationships between DNAm across diverse genomic contexts and gene regulation (Jones, 2012). In the present study, sex-associated DMCs were significantly depleted in promoters and enriched in intergenic regions of the genome. It is possible that intergenic differential methylation patterns reflect distal cisregulatory interactions as has been demonstrated in mammalian systems (Jones, 2012;Schmidl et al., 2009;Stadler et al., 2011;Stroud et al., 2011;Wiench et al., 2011). For example, human immune cells exhibit cell-type-specific DNAm patterns in genomic regions distal to promoters that are associated with methylation-dependent enhancer activity (Schmidl et al., 2009). Such DNAm patterns also correlate with histone modifications including histone H3 lysine 4 methylation (Schmidl et al., 2009). It is intriguing to consider how During TSD, the gonadal DNA methylome undergoes widespread remodelling in response to incubation temperature. These epigenomic changes probably play roles in both the establishment of sexual fate during the discrete window of development in which sex determination is sensitive to temperature as well as the maintenance of sexual fate into adulthood (Navarro-Martín et al., 2011;Weber & Capel, 2021). Compared to temperature-associated DNAm patterns in embryonic gonads immediately following sex determination, temperature influences on the DNA methylomes of hatchling blood cells were considerably less extensive. Further, the loci acquiring temperature-associated DNAm patterns appeared largely unique to each tissue. This is perhaps not surprising considering DNAm patterns play a key role in the maintenance of cellular identity (Bock et al., 2012;Ziller et al., 2013) and thus tend to diverge markedly depending on tissue and cell type (Blake et al., 2020;Christensen et al., 2009).
Nonetheless, previous studies have identified subsets of loci displaying consistent DNAm responses to environmental cues across tissues, termed "metastable epialleles" (Bourc'his, 2021;Rakyan et al., 2002). In a study of European sea bass, nine genes showed consistent differential methylation related to past developmental temperature across all tissues examined, including brain, muscle, liver and testis (Anastasiadi et al., 2021). While results from the present study suggest DNAm responses to incubation temperature are largely tissue-specific, it should be noted that our  Barua et al., 2014;Borghol et al., 2012;Dolinoy et al., 2007;Heijmans et al., 2008). Even within the embryonic period, certain developmental windows are characterized by unique environmental sensitivity (Faulk & Dolinoy, 2011). In the context of TSD, the thermal plasticity of gonadal fates is limited to a specific period of development termed the thermosensitive period during which sexually dimorphic epigenomic patterns are established in response to temperature (Ge et al., 2018;McCoy et al., 2015). These gonadal DNAm patterns are then maintained across mitotic cell divisions. However, somatic tissues are likely to differ in their pattern of epigenetic thermal plasticity based on multiple factors including cell type-specific transcriptional activity F I G U R E 7 Comparison of temperature-associated DNA methylation patterns in blood vs. embryonic gonad just after temperaturedependent sex determination. (a) Proportion of covered CpG sites shared among blood and gonad RRBS data (n = 181,834) that are differentially methylated with respect to incubation temperature. The p-value is derived from a Fisher's exact test to determine if the proportion of differentially methylated vs. nondifferentially methylated loci is different between tissues. (b) Comparison of per cent methylation difference (i.e., effect size) of temperature-associated differentially methylated cytosines (DMCs) identified in each tissue. (c) Comparison of the distribution of temperature-associated DMCs identified in each tissue with respect to gene context. Asterisk indicates a significant difference in proportion of DMCs in a context category compared to the shared covered cytosines at an FDR threshold of 10%. The only significant differences we observed were depletion of gonadal temperature-associated DMCs in promoter regions and enrichment of gonadal temperature-associated DMCs in introns.  (Almstrup et al., 2016;Hanna et al., 2012;Nugent et al., 2015;Zhao et al., 2017). Interestingly, many sex-associated DMCs identified in this study were unique to individuals from the 34.5°C treatment. It is possible that in the absence of incubation temperature-associated DNAm variation, sex-associated DNAm patterns are more readily detected. Alternatively, this finding raises the possibility that females and males from incubation temperatures yielding mixed sex ratios may exhibit a greater degree of sexual dimorphism at the level of the DNA methylome compared to those individuals from incubation temperatures that only yield one sex. Given that this observation stems from a single incubation temperature and relatively small sample size, further work examining patterns of interand intrasexual variation in DNAm in individuals from temperatures producing mixed sex ratios is needed to test the degree to which epigenetic sexual dimorphism varies with the developmental thermal environment.
Overall, sex-associated DNAm patterns in blood cells represent promising candidates to serve as nonlethal indicators of phenotypic sex in early-life stages of TSD species. However, results regarding the predictive power of the DNAm-based model presented here should be interpreted with caution given the limited sample size. Further, genome-wide approaches to profiling DNAm are not currently a practical option for conservation efforts to monitor hatchling sex ratios in the field. The utility of these DNAm marks as a conservation tool depends on validation using more targeted approaches (e.g., bisulphite sequencing PCR) across many samples originating from multiple geographical sites. Given the links between genomic sequence variation and epigenomic variation (Fargeot et al., 2021;Foust et al., 2016;Wogan et al., 2020), incorporating samples from distinct populations is crucial to develop a generalizable tool for use in conservation contexts. Further, a quantitative comparison of the reliability and cost-efficiency of different potential molecular markers of sex, for example blood DNAm patterns, sex-biased protein levels (Tezak et al., 2020), and high-sensitivity steroid hormone quantitation (Bae et al., 2021), in TSD species is warranted.
In the face of rapidly changing environmental conditions, monitoring the responses of key demographic parameters, such as offspring sex ratios, to climate change in real time will play a critical role in the conservation of thermally sensitive species. This is the first study to use somatic DNAm patterns to predict phenotypic sex in a species with environmental sex determination. Intriguingly, we also demonstrate that DNAm patterns hold the potential to provide information about past incubation temperatures, though this finding requires validation in the context of fluctuating incubation temperatures as well as in older animals. Understanding how genetic and environmental variation contribute to these methylation patterns represents a critical target for future research that will advance efforts to realize the full potential of epigenomic tools in conservation contexts.

CO N FLI C T O F I NTE R E S T
The authors declare no competing interests.

DATA AVA I L A B I L I T Y S TAT E M E N T
The raw sequencing reads and associated metadata from reducedrepresentation bisulphite sequencing libraries are available via the NCBI Sequence Read Archive (SRA) under the BioProject PRJNA857046 (whole blood) and PRJNA857084 (embryonic gonad).

D I SCL A I M ER
This report was prepared as an account of work sponsored by an agency of the United States Government. Neither the United States Government nor any agency thereof, nor any of their employees, makes any warranty, express or implied, or assumes any legal liability or responsibility for the accuracy, completeness, or usefulness of any information, apparatus, product, or process disclosed, or represents that its use would not infringe privately owned rights. Reference herein to any specific commercial product, process, or service by trade name, trademark, manufacturer, or otherwise does not necessarily constitute or imply its endorsement, recommendation, or favouring by the United States Government or any agency thereof.
The views and opinions of the authors expressed herein do not necessarily state or reflect those of the United States Government or any agency thereof.