Authors: Jessica M. Judson (1Department of Ecology, Evolution, and Organismal Biology, Iowa State University, Ames, IA 50011, USA; 3Current Address: W. K. Kellogg Biological Station, Departments of Fisheries and Wildlife & Integrative Biology, Michigan State University, Hickory Corners, MI 49060, USA), Luke A. Hoekstra (1Department of Ecology, Evolution, and Organismal Biology, Iowa State University, Ames, IA 50011, USA; 2Current Address: Department of Integrative Biology, Oklahoma State University, Stillwater, OK 74078, USA), Fredric J. Janzen (1Department of Ecology, Evolution, and Organismal Biology, Iowa State University, Ames, IA 50011, USA; 3Current Address: W. K. Kellogg Biological Station, Departments of Fisheries and Wildlife & Integrative Biology, Michigan State University, Hickory Corners, MI 49060, USA)
Categories: Article, Range expansion, surfing, painted turtle, ectotherm, population genomics
Source: Molecular ecology
Doi: 10.1111/mec.17269
Authors: Jessica M. Judson, Luke A. Hoekstra, Fredric J. Janzen
Environmental conditions vary greatly across large geographic ranges, and yet certain species inhabit entire continents. In such species, genomic sequencing can inform our understanding of colonization history and the impact of selection on the genome as populations experience diverse local environments. As ectothermic vertebrates are among the most vulnerable to environmental change, it is critical to understand the contributions of local adaptation to population survival. Widespread ectotherms offer an opportunity to explore how species can successfully inhabit such differing environments and how future climatic shifts will impact species’ survival. In this study, we investigated the widespread painted turtle (Chrysemys picta) to assess population genomic structure, demographic history, and genomic signatures of selection in the western extent of the range. We found support for a substantial role of serial founder effects in shaping population genomic demographic analysis and runs of homozygosity were consistent with bottlenecks of increasing severity from eastern to western populations during and following the Last Glacial Maximum, and edge populations were more strongly diverged and had less genetic diversity than those from the center of the range. We also detected outlier loci, but allelic patterns in many loci could be explained by either genetic surfing or selection. While range expansion complicates the identification of loci under selection, we provide candidates for future study of local adaptation in a long-lived, widespread ectotherm that faces an uncertain future as the global climate continues to rapidly change.
How do certain species inhabit entire continents, despite the extreme differences in environmental conditions experienced by populations across space? To tolerate these considerable shifts in environment, many populations have locally adapted, such that transplanting individuals from one population to another would result in a decrease in their fitness compared to native individuals (Hereford, 2009). Particularly when migration is uncommon, populations are expected to locally adapt toward a fitness optimum for their environment, independent of the selective pressures faced by other populations (Orr, 1998). Examples of local adaptation are numerous (see reviews by Fraser et al., 2011; Leimu & Fischer, 2008; Sanford & Kelly, 2011; Urban et al., 2014), and local adaptation is often invoked to explain persistence across large geographic ranges.
Individuals can also respond to diverse environmental challenges through phenotypic plasticity (Bodensteiner et al., 2023; Bonamour et al., 2019; Torres-Dowdall et al., 2012). In long-lived iteroparous species, phenotypic plasticity is an important mechanism for adjusting to diverse environmental conditions experienced within an individual’s lifetime (Gunderson & Stillman, 2015; Nussey et al., 2007; but see Radchuk et al., 2019). However, in the face of anthropogenic climate change, the capacity of species with long generation times for rapid adaptive evolution is less certain (Moritz & Agudo, 2013). To understand both the historical evolutionary pressures that shaped species’ current ranges and the future risk of extinction from anthropogenic climate change, it is essential to assess the contributions of both phenotypic plasticity and local adaptation to phenotypic change across populations of long-lived species (Valladares et al., 2014). New genomic resources allow the investigation of patterns of selection across the genome, past demographic patterns and current population genetic structure (Savolainen et al., 2013). This new information can contribute to our understanding of how long-lived species with widespread geographic distributions inhabit such a wide range of natural environmental conditions and how they may respond to future climate change (Meek et al., 2023).
The predicted survival outcome for many reptiles appears bleak due to the rapid pace of climatic shifts exacerbated by habitat loss and decreased population sizes (Böhm et al., 2016; Gibbons et al., 2000; Sinervo et al., 2010). Thermosensitive traits such as temperature-dependent sex determination (TSD) and anoxia tolerance make these species particularly vulnerable. Temperature changes of only a few degrees could destabilize sex ratios, especially in small populations (Mitchell & Janzen, 2010; Mitchell et al., 2010), or influence individuals’ ability to prepare adequately for warming winter conditions (Moss & MacLeod, 2022). In reptiles, most studies investigating potential responses to climate change suggest that plasticity plays an important, and perhaps predominant, role in mediating organismal responses to temperature (e.g., Janzen et al., 2018; Refsnider & Janzen, 2012; Refsnider & Janzen, 2016; Urban et al., 2014). Although plasticity is likely important for the survival of reptiles with long generation times, the potential contribution of local adaptation to future survival of reptile populations is perhaps less understood than the contribution of phenotypic plasticity (Urban et al., 2014). We currently have a limited grasp of the amount of standing genomic variation in long-lived reptile populations, how that variation has been shaped by past demographic events (e.g., population expansions, bottlenecks), and the degree to which this variation might contribute to adaptation. Investigating population genomic structure, demographic history, and genomic evidence of selection in long-lived reptiles can inform management of these charismatic species and help predict future responses to climate change (Macdonald et al., 2018).
Painted turtles (Chrysemys picta) – long-lived ectotherms with TSD - are an excellent reptile to evaluate genomic variation and signatures of selection across a large geographic range. Painted turtles are widespread across North America, with populations subjected to disparate macro-environmental conditions. For example, in the species’ western extent, maximum air temperatures experienced among localities range from <26°C to 35°C and minimum air temperatures range from 0°C to <−22°C. Divergence in color, scute patterns, nesting behavior, adult body size, reproductive effort, characteristics of TSD reaction norms (e.g., transitional range of temperatures), and thermal reaction norms for hatchling body mass and incubation time have been documented across populations (Bodensteiner et al., 2019; Bodensteiner et al., 2023; Carter et al., 2019; Edge et al., 2017; Ernst & Lovich, 2009; Iverson & Smith, 1993; Janzen et al., 2018; Ultsch et al., 2001). Phenotypic plasticity explains some variation in phenotypes among populations (e.g., date of first nesting, Janzen et al., 2018; nesting behavior, Refsnider & Janzen, 2012). However, population genomic structure and genomic signatures of selection have not yet been assessed in this species. In this study, we investigated population genomic structure, demographic history, and genomic signatures of selection in seven locations spanning the western range of painted turtles. We also assessed patterns of divergence in particular genes associated with TSD and extreme anoxia tolerance (e.g., Fanter et al., 2020; Ge et al., 2018), as we hypothesized these two traits may be under differential selection across painted turtle populations experiencing different temperature extremes. Understanding both demographic history and selection dynamics in a widespread, long-lived ectotherm provides insight into how species successfully colonize vast geographic ranges.
Painted turtles (Chrysemys picta) have a unique distribution among turtles, with the widest natural range of any North American turtle species (Ernst & Lovich, 2009). The current distribution encompasses the entire eastern coast of the United States and stretches west across the northcentral United States and southern Canada, with a patchy distribution in the southwestern United States (Fig. S1; Bodensteiner et al., 2019; Carter et al., 2019). Fossil evidence suggests the species predates the Pleistocene glaciation period, and limited genetic evidence points to multiple range contractions as glaciers expanded (Starkey et al., 2003). Using one mitochondrial locus, Starkey et al. (2003) found that the southeastern populations were more divergent from other populations and were likely the source of multiple colonizations of the northern United States and Canada following glacial retreat. Further, very low genetic divergence was detected across northern and southwestern populations, particularly in the western half of the species’ range, despite the species’ age and widespread distribution. This result suggests that glaciation forced a range contraction in the east, followed by a rapid westward radiation of painted turtles. Additionally, the reference genome for painted turtles shows slow rates of sequence evolution compared to other vertebrate lineages, which may further compound a putative lack of genetic diversity observed among populations (Shaffer et al., 2013). More recent studies of genetic structure in painted turtles, again with limited number of markers, have consistently reported low divergence across much of the range (Jensen et al., 2015), including a few western populations (Reid et al., 2019). The study by Reid et al. (2019) corroborates the hypothesis that the western populations are the result of a rapid radiation westward following glaciation. Taken together, previous work on the genetic structure of this species has relied on a small number of neutral markers and suggests climate has strongly influenced the past and current range of painted turtles.
Though low genetic diversity has been consistently reported across western painted turtle populations, phenotypes do vary across locations from the species’ western range. For example, body size and age at maturity of adult painted turtles was positively correlated with latitude and negatively correlated with mean annual temperature (Iverson & Smith, 1993). Additionally, thermal reaction norms for painted turtle hatchling phenotypes varied in common-garden conditions; responses of incubation duration and hatchling mass to temperature varied among hatchlings from different locations, but in a pattern inconsistent with latitude (Bodensteiner et al., 2019). Temperature and precipitation patterns vary widely across the species’ range in the western United States (Table S1), which likely exerts different selective pressures on these ectotherms, particularly at thermal extremes. As a model of a long-lived, widespread ectothermic amniote with an exceptional fossil record (see Starkey et al., 2003), a genome-wide study of genetic differentiation in painted turtles is warranted to investigate demographic patterns and putatively adaptive variation.
We collected clutches of painted turtle hatchlings or eggs from seven locations across the western 2/3 of the species’ range (Table 1; Fig. S1) and incubated eggs at Iowa State University (ISU) until hatching (see Bodensteiner et al., 2019; Refsnider et al., 2014). We euthanized hatchlings with a pericardial overdose of 0.5 mL of 1 sodium water, assigned sex according to macroscopic examination of gonads (Schwarzkopf & Brooks, 1985; Warner et al., 2014), and stored hatchlings in 70% ethanol until DNA extraction. All research followed ISU IACUC protocols (6-08-6583-J and 12-03-5570-J) with appropriate permits from the associated agencies of each sampling location (see Data Accessibility). To avoid any bias of relatedness, we genotyped only one hatchling from each clutch, with sample sizes ranging from 5 to 38 hatchlings from each location (Table 1). We extracted DNA from hatchling liver tissue using a DNeasy Blood and Tissue Kit (Qiagen, Germantown, MD) or a phenol-chloroform DNA extraction protocol (Sambrook et al., 1989; see SI). We quantified DNA with a NanoDrop^™^ 2000 Spectrophotometer and assessed DNA purity and quality using a 1% agarose gel before we performed DNA sequencing on 164 individuals.
We used a restriction site-associated DNA sequencing (RADseq) approach with a single-restriction enzyme, MseI, ligated resulting fragments with two adapters (6bp or 8bp barcode), and pooled samples for size selection. We selected fragments of 275-300bp with gel electrophoresis and confirmed size with an Agilent^®^ 2100 BioAnalyzer before paired-end 150bp sequencing on an Illumina HiSeq platform (Illumina, San Diego, CA). Following sequencing, we removed raw reads that contained more than 10% unknown bases, removed reads with more than 50% low quality bases (Q ≤ 5), and trimmed adapter sequences with cutadapt v2.5 (Martin, 2011). Finally, we used the BWA-mem algorithm (Li, 2013) in bwa v0.7.17 (Li & Durbin, 2009; Li & Durbin, 2010) to align reads to the painted turtle reference genome v3.0.3 downloaded from NCBI, which was generated with a female painted turtle from southern Washington state, USA (Shaffer et al., 2013).
We processed aligned reads for variant calling by removing unmapped and improperly paired reads, reads not in primary alignment, and alignments with mapping quality (MAPQ) < 30 using SAMtools v1.9 (Li et al., 2009). To call variants, we used a modified Sentieon DNAseq workflow (v201808.01; https://support.sentieon.com/manual/DNAseq_usage/dnaseq/). We used Sentieon’s Haplotyper algorithm to generate genomic variant call format files (GVCFs), including only sites with a minimum base quality > 20, and generated a VCF with raw variant calls across all samples with the GVCFtyper algorithm. We created a mask file containing both confident variant and invariant reference sites using ‘–emit_mode confident’. Following genotyping, we applied hard filters to single nucleotide polymorphisms (SNPs) and indels/mixed sites separately using GATK’s VariantFiltration and SelectVariants (v4.0.4.0, DePristo et al., 2011; McKenna et al., 2010). We removed SNPs if they fit any of the following low variant confidence (QD < 2.0), strand bias (FS > 60.0 or SOR > 4.0), low mapping quality (MQ < 40.0), large differences between reference and variant mapping quality (MQRankSum < −12.5), or bias in position of alleles within reads (ReadPosRankSum < −8.0; e.g., Pfeifer et al., 2018). We removed indels if QD < 2, FS > 200, or ReadPosRankSum < −20.
We further filtered the variant and invariant VCFs using VCFtools v0.1.14 (Danecek et al., 2011). We removed individuals with low average depth (< 2X), excluded variants with mean genotype quality < 20.0, and restricted to biallelic sites where applicable. We then removed indels (filtered above), as well as 3bp upstream and downstream of the indel, from both the SNP and the invariant sites VCFs. Next, we removed SNPs that deviated strongly from Hardy-Weinberg equilibrium (HWE; P < 10^−7^). As the sampled locations are extremely geographically distant, we set this threshold low to reduce the likelihood of removing real variants generated by substantial population genomic structure (e.g., Martin et al., 2018). We removed sites fixed for the non-reference allele and both variant and invariant sites with >25% of genotypes missing across individuals. This procedure yielded 41,546,386 genotyped sites across the painted turtle genome (genome size ~2.6Gb; Shaffer et al., 2013). We then applied a stricter set of hard filters for population genomic analyses (Pfeifer et al., 2018), which are outlined here. First, we assigned SNP genotypes with a depth less than 4 or greater than 100 as missing to avoid including any low coverage reads or paralogous regions. Second, we partitioned the painted turtle genome scaffolds into 1.5kb blocks and assigned all SNPs and monomorphic sites to their corresponding blocks using a custom R script (v3.6.3; R Core Team, 2020) with “dplyr” v0.7.5 (Wickham et al., 2018). We removed any blocks with a median distance between SNPs less than 3 or with less than 150 genotyped sites (both variant and invariant) to prevent the inclusion of SNPs in repetitive or misaligned regions. This approach produced a dataset with 37,832,754 invariant and 566,852 variant sites over 153,416 blocks, which was used for calculating nucleotide diversity and percentage of polymorphic loci.
For analyses using only variant sites, we further reduced missing data among SNPs and individuals by removing SNPs with greater than 10% missing genotypes and confirmed that individuals had greater than 75% SNPs genotyped in the resulting dataset. The resulting VCF file contained 253,365 SNPs across 161 individuals. Then, we applied a minor allele frequency (MAF) filter such that singletons were excluded, and we retained only one SNP with the least missing data per 1.5kb block to reduce linkage among markers, for a final dataset of 48,060 SNPs used for population genomic analyses. A summary of filtering steps and number of sites can be found in Table S2, and all custom scripts can be found in the repository (Judson et al., 2023).
To calculate intrapopulation (π) average nucleotide diversity values, we used VCFtools v0.1.14 --site-pi (Danecek et al., 2011) to calculate per-site nucleotide diversity for each population and then weight-averaged the estimate by the number of genotyped sites following filtering (38,399,606; Martin et al., 2018). We used the same set of sites and weight-averaging to calculate interpopulation (dxy) nucleotide diversity (with R package “PopGenome” v2.7.5; Pfeifer et al., 2014) and percentage of polymorphic loci across sampling locations. We used the final SNP dataset including all filters (48,060 SNPs) to calculate summary statistics for each sampling location using “vcfR” v1.12.0 (Knaus & Grünwald, 2017), “poppr” v2.8.7 (Kamvar et al., 2014), “adegenet” v2.1.3 (Jombart & Ahmed, 2011), “heirfstat” v0.5.7 (Goudet & Jombart, 2020) and “SNPRelate” v1.22.0 (Zheng et al., 2012) in R. These statistics include average percentage of missing genotypes per individual, number of private variable sites, expected (HE) and observed (HO) heterozygosity, average inbreeding coefficient (FIS), and Weir and Cockerham’s overall and pairwise FST (Weir & Cockerham, 1984).
To understand population genomic variation, we performed a principal components analysis (PCA) in R using “SNPRelate” v1.22.0 (Zheng et al., 2012) downloaded from BiocManager v1.30.10 (Morgan, 2019). We further assessed groupings in our data using the program STRUCTURE v2.3.4 (Pritchard et al., 2000) with the admixture model and the following 100,000 burn-in and 100,000 repetitions for each run of K from K=1 to K=10, correlated allele frequencies among populations, the alternative, population-specific ancestry prior (which entails allowing a separate alpha for each population; POPALPHAS=1), and an initial alpha value of 0.14, where alpha is the result of 1/K and K is the assumed number of populations. These settings were chosen to account for unequal sampling sizes from the seven locations (Wang, 2017). We chose K=10 to allow some subdividing within sample locations (N=7), but did not expect population structure within sampling locations given the focused spatial sampling. Each value of K was run 20 times using the above parameters, as recommended by Gilbert et al. (2012). To select the optimal K from the STRUCTURE runs, we used multiple methods (Janes et al., 2017): we assessed the posterior probabilities for each K (Pritchard et al., 2000) and plotted Ln Pr(X|K) (Pritchard & Wen, 2003) and ΔK (Evanno et al., 2005) using CLUMPAK v. 1.1 (Kopelman et al., 2015) and STRUCTURE HARVESTER (Earl & Vonholdt, 2012).
We built a maximum likelihood population tree using TreeMix v1.13, which uses a graph-based approach to model both gene flow and population splits (Pickrell & Pritchard, 2012). We used 500 bootstrap replicates, and assumed no migration given the distance among sampled locations. We tested for isolation by distance (IBD) using a Mantel test comparing straight-line geographic distance (km) between locations to pairwise genetic differentiation, FST/(1-FST) (Rousset, 1997), implemented in “vegan” v.2.5.7 (Oksanen et al., 2020) with 10,000 permutations. We also assessed correlations between pairwise distance and differences in annual mean temperature (BIO1) and precipitation (BIO12; Table S1) between locations using a Mantel test, as strong correlations could influence our attribution of population structure to demographic history and selection.
Runs of homozygosity (ROH) can be informative of both past demographic changes (e.g., bottlenecks) and recent inbreeding within populations (reviewed in Ceballos et al., 2018). Specifically, long tracts of homozygosity are generated by recent inbreeding among individuals, while more numerous short ROH can be indicative of population bottlenecks farther back in time. We assessed runs of homozygosity in the variant dataset with all filters except the final filter retaining one SNP per block using PLINK v1.9 (Purcell et al., 2007). Before analysis, we separated individuals by sampling location and removed SNPs that were fixed within each location as recommended in Shafer et al. (2016) using GATK (v4.0.4.0, DePristo et al., 2011; McKenna et al., 2010). We used the –homozyg algorithm in PLINK with the following window size of 25 SNPs, defined windows as homozygous if one or less heterozygous sites were found in the window, allowed no more than 5 missing sites within a window, and set the window threshold for a SNP being defined as in a homozygous segment of a chromosome (hit rate) at 0.05. For a region to be defined as a ROH, segments must include ≥ 25 SNPs, cover ≥ 1,000 kb, contain a maximum of one heterozygous site, contain minimum SNP density of 1 per 50 kb, and have maximum distance between two adjacent SNPs of ≤ 1,000 kb (Grossen et al., 2018).
Finally, we used StairwayPlot v2.1.1 (Liu & Fu, 2015; Liu & Fu, 2020) to infer the demographic history of painted turtles from the four sampled locations with a larger sample size (N ≥ 27). We used the variant dataset without the MAF filter or the filter retaining one SNP per block (253,365 SNPs). To create an equivalently filtered monomorphic site comparison for correctly assessing the total number of genotyped sites for demographic analysis, we further filtered the invariant site VCF to remove sites with greater than 10% missing genotypes (Table S2). We created population-specific VCFs using VCFtools v0.1.14 (Danecek et al., 2011). To create population-specific folded site-frequency spectra (SFS), we used the python script easySFS (Gutenkunst et al., 2009; Overcast, 2023). We downsampled to maximize the number of SNPs retained and set --total-length to the sum of variant and invariant sites (26,707,837). We used a generation time of 11 years (Reid et al., 2019) and a mutation rate per site per generation of 4.61 x 10^−9^. This painted turtle mutation rate was based on parent-offspring trio sequencing of turtles from the Illinois sampling site in this study (Bergeron et al., 2023).
To detect outlier loci potentially linked to selection across the sampled locations, we first applied the Bayesian approach of BayeScan v2.1 (Foll & Gaggiotti, 2008). For this analysis, we used the variant dataset that contained all filtering steps described above except the final filter retaining one SNP per block (212,670 SNPs; Table S2). We included these potentially linked SNPs to give increased support for any detected outlier SNPs; if selection is causing a deviation from the population mean FST, the surrounding SNPs should deviate in a similar manner (Storfer et al., 2018). We used a burn-in of 50,000, 100,000 iterations, and prior odds of 100 for the neutral model to reduce the number of false positives. After running BayeScan twice, we evaluated convergence using the R package CODA (Plummer et al., 2006) and identified any loci with a q-value less than 0.05. We performed two additional BayeScan runs with the same settings on a dataset excluding Kansas turtles and two admixed Oregon turtles (hereafter, “reduced dataset”; see Results), given the small sample size from Kansas.
We also used a separate approach to detect outlier loci, “pcadapt” v4.3.3 (Privé et al., 2020), which uses a PCA approach to find markers that are significantly related to population structure. After confirming that linkage disequilibrium did not influence outlier detection (Luu et al., 2020), we used the same dataset as the BayeScan analysis including 212,670 SNPs. After assessing values of K from 1–20, we used Cattell’s rule to select K=5; higher PCs split individuals from the same sampling location (Fig. S2). Outliers were determined using two p-value cutoff q-values with α = 0.05 (similar to BayeScan) and using Bonferroni correction with α = 0.001. We ran the same analyses with the reduced dataset. While we set the threshold for the HWE filter to be low (P < 10^−7^), the use of an HWE filter can exclude loci that deviate due to population structure and that may be under selection. Thus, we added loci previously excluded for excess homozygosity (see SI) to previously included loci (217,821 total), excluded Kansas and admixed Oregon turtles, and used pcadapt to detect outliers (Bonferroni correction α = 0.001). We assessed the genomic location, including whether outliers were in genes or exons, for outlier SNPs using BEDTools v2.27.1 (Quinlan & Hall, 2010). We used VCFtools v0.1.14 (Danecek et al., 2011) to calculate FST for all loci included in outlier detection to assess patterns across the genome. Finally, we assessed whether the gene ontology (GO) functional category of genes represented in the outlier list generated from pcadapt were statistically overrepresented when compared to all represented genes in the SNP dataset using PANTHER v17.0 with FDR correction (Mi et al., 2019; Thomas et al., 2022).
Given painted turtles exhibit two remarkable aspects of biology, TSD and extreme anoxia tolerance (Janzen, 1994; Ultsch & Jackson, 1982), that conceivably could be linked to adaptation to macro-environmental conditions, we further explored potential adaptation at the candidate gene level. First, we queried coverage of genes in our RADseq dataset that are implicated in the sex-determination mechanism, as these genes may play a role in differentiation of the transitional range of temperatures across these populations (Carter et al., 2019). The genes we chose to assess for variation were Doublesex and Mab-3–Related Transcription factor 1 (DMRT1), the expression of which was recently confirmed to be the key sex-determining factor in birds (Ioannidis et al., 2021), and lysine demethylase 6B (KDM6B), a gene that promotes transcription of DMRT1 in red-eared slider turtles (Ge et al., 2018). Second, we assessed variation in genes identified as differentially expressed during anoxia or re-oxygenation in painted turtles (Fanter et al., 2020). We plotted results of all analyses using “ggplot2” v3.3.3 (Wickham, 2016) in R.
In total, 164 painted turtle hatchlings representing 7 sampling locations were sequenced for this project (Table 1; NCBI BioProject PRJNA789046, Judson et al., 2023), with an average of 3.77 million high-quality reads per individual and 3.5 million reads mapped to the reference genome. Joint genotyping yielded 3.53 million SNP loci; individual depth filters and site filters resulted in a final dataset of 48,060 SNPs from 161 individuals with an average of 4.43% missing genotypes across individuals (Table 1, Table 2, Table S2). With such a large number of SNPs, a representation of genetic diversity to understand population genetic structure in these locations was likely captured even in locations where the number of individuals sampled was small (Nazareno et al., 2017).
Summary statistics suggested that genetic variation varied widely among populations. Polymorphism ranged from 0.48% to 1.23% (Table 2), with the lowest values of polymorphism and private alleles in the northwestern populations and the highest in the easternmost (Illinois) site. The Illinois population had much greater genetic variation than other sampled locations, with private alleles an order of magnitude more numerous than any other population. Weir and Cockerham’s overall weighted FST across populations was 0.136 and average FST across loci was 0.089.
Pairwise FST, dxy, PCA, and STRUCTURE results all depicted a similar pattern of genetic structure across the sampled populations. Pairwise FST values suggested moderate to high differentiation among many populations (Table 3). However, Illinois, Nebraska, Minnesota, and Kansas had lower levels of differentiation than other pairwise comparisons, and Oregon and Idaho pairwise FST was similarly low. The first three principal components (PC) explained 20% of the variation in the dataset, with percentages of 9.4%, 6.2% and 4.5% of the variance explained, respectively. The first PC was consistent with geographic positioning of populations from east to west, with Illinois and the northwestern individuals farthest from each other (Fig. 1A). The second PC separated the New Mexico individuals from the more northern populations, and the third PC separated New Mexico individuals from Nebraska, Kansas, and Minnesota individuals (Fig. S3). The results of the various estimators of best K for the STRUCTURE runs suggested six genetically distinct groups were sampled (Fig. 1C; Fig. S4–S8). All sampled locations were distinct in the STRUCTURE plot except individuals from Oregon and Idaho, which grouped as a single population. ΔK suggested the optimal K was 2 (Fig. S9), but this result is not convincing given the frequent recovery of K=2 using this method (Janes et al., 2017) and other estimators suggesting K=6. Painted turtles from the Oregon and Idaho locations were similar in allele frequencies; these individuals overlapped in the PCA and displayed low pairwise FST and dxy (Table 3). Two individuals from Oregon were distinguished from the other Oregon individuals in both the PCA and STRUCTURE analyses and may be the result of admixture between Oregon and midwestern populations (IL, MN; Fig. 1). We calculated individual relatedness using PLINK, and both individuals appeared not to be closely related (see SI).
The population tree from TreeMix analysis explained 99.8% of the covariance in allele frequencies among locations. There was substantial drift among populations, with Oregon and Idaho grouping closely and the remaining populations showing branching consistent with their east-west location (Fig. 2). The Mantel test of IBD supported a significant relationship between pairwise geographic distance and pairwise genetic distance (P = 0.012; Fig. 1B). Mantel tests indicated no correlations between pairwise distance and differences in annual mean temperature (P = 0.66) and precipitation (P = 0.89).
We assessed the demographic history of western painted turtles using both patterns of ROH and the site-frequency spectrum. We found that the total length of ROH was correlated with the number of ROH (R^2^=0.9763; Fig. 3A). The relationship was linear, with individuals from the most genetically variable population (Illinois) having fewer ROH, individuals from Oregon and Idaho having more of their genome in ROH and more numerous ROH, and other populations displaying intermediate ROH numbers consistent with their geographic distance from the Illinois population (Fig. 3A). ROH were then partitioned into either small and intermediate ROH < 5 megabases (MB) or long ROH > 5 MB, the latter of which indicate recent inbreeding (e.g., Grossen et al., 2018; McQuillan et al., 2008). When comparing long ROH, populations did not vary in the total length of long ROH in the genome, as would be expected under a scenario of increased inbreeding among individuals within certain populations (Fig. 3B). Instead, populations varied in the total length of the genome that resides in short and intermediate ROH, with Oregon and Idaho individuals having much more numerous ROH < 5MB than other populations (Fig. 3C). This result is more consistent with past bottlenecks having generated the observed ROH than with recent inbreeding, particularly when coupled with the linear increase of both ROH length and number with geographic distance among populations (Fig. 3A).
To further investigate whether bottlenecks could explain the patterns of ROH across populations, we used StairwayPlot2 to assess demographic history in the four sampling locations with larger sample sizes (Illinois, Minnesota, Nebraska, and Idaho, Table 1). We found that the western populations in Nebraska and Idaho show clear evidence of a bottleneck occurring during and after the Last Glacial Maximum (LGM; Fig. 4, Fig. S10, Fig. S11), consistent with the ROH results above. Additionally, the severity of the bottleneck in the Idaho population is much greater than that in Nebraska, consistent with serial founder effects during range expansion. Effective population size estimates are large, but the average values match those of the combined Chrysemys picta bellii group in Reid et al. (2019). While the mutation rate is directly from a population used in this study (Bergeron et al., 2023), generation time could vary in painted turtle populations (Spencer & Janzen, 2010; Wilbur, 1975), and we cannot account for variation in sex ratios, so timing of bottlenecks and estimates of effective population size may not be precise.
After validating convergence, BayeScan identified 43 FST outlier loci at the 0.05 FDR for the dataset including all individuals (Fig. S12, Table S3) and 53 outlier loci for the reduced dataset (Fig. S13, Table S4; discussed here). When we investigated the location of these SNPs in the genome for proximity to genic regions, 34 of the 53 outliers corresponded to genic regions, and two of those were located in a coding region according to the BEDTools query. According to Ensembl’s Variant Effect Predictor (VEP; release 110, Cunningham et al., 2022), one variant was downstream of the LDLRAD3 gene, and the other variant was a synonymous SNP in the RBBP6 gene. We investigated genotypic patterns among populations at the 53 outlier loci and found that all 53 outlier loci were segregating in the Illinois population. Further, these loci were often fixed in the sampling sites farthest away from Illinois (Oregon, Idaho, and New Mexico populations).
The second outlier detection method (pcadapt) with q-value α = 0.05 (similar to BayeScan) detected 2,802 outliers with the full dataset and 2,597 outliers with the reduced dataset. These outliers overlapped with the BayeScan results at ten SNPs using the full dataset (Table S5) and 13 in the reduced dataset (Fig. S14, Table S6), though some of the overlapping SNPs were different. Genes that overlapped in both outlier lists with BayeScan outliers were TNR, ZNF436 – like, CASKIN, and SIM2. We focused on a narrower set of outlier loci using Bonferroni correction with α = 0.001 and found 228 outliers in the reduced dataset (Fig. S14). Plotting FST across the genome (Fig. 5, Fig. S14) did not show a concentration of high FST or outliers in any specific genomic regions. Outliers detected with pcadapt tended to include the SNPs with the highest FST, while BayeScan outliers had lower FST. Including loci previously excluded with the HWE filter resulted in 476 outliers. New outliers from this SNP set had the highest FST of all loci (Fig. 5), confirming that the HWE filter removed loci strongly differentiated among populations. 228 of these SNPs (48%) were in 166 unique genes, and 10 were within exons. Of the SNPs that were not downstream of a gene according to VEP, four were synonymous mutations and three were missense mutations (in ZAR1l, AKAP6, and SCML2). The missense mutation in ZAR1l appears to be the result of a de novo mutation in the New Mexico population with a high frequency (0.75). Genotypic patterns were more varied than in the BayeScan some loci were fixed in the Illinois population and variable in populations farther west. The overrepresentation test in PANTHER found no statistically significant GO category overrepresented in the outlier SNP genes (Fig. S15, Fig. S16).
We also assessed variation at specific genes implicated in TSD and anoxia tolerance in painted turtles. Although we did not find coverage of KDM6B, DMRT1 was sequenced in our dataset with 22 SNPs. None of these SNPs were in coding regions, and these SNPs did not appear to be correlated with sex of sampled hatchlings, though most hatchlings were incubated in field conditions and thus we cannot assess the influence of genotype on sex at the pivotal temperature of ~28°C (when theoretical sex ratio is 1; Telemeco et al., 2013). Second, we had coverage of ANKRD1, ATF3, CA1, SYR61, ENDOD, PTGS2, S100A1, and TNC, which are all identified as differentially expressed during anoxia or re-oxygenation in painted turtles (Fanter et al., 2020). One SNP in CA1 was located in a coding region, and changed an amino acid from proline to glutamine. This SNP was found in five Illinois, seven Minnesota, and one Nebraska individual. Four other SNPs associated with amino acid changes in coding regions of CYR61 and TNC were found at low frequency in one to three individuals.
With genomic data, we can now investigate in greater detail the genomic architecture of local adaptation among populations of widespread species and understand past demographic changes that led to species distributions. In this study, we used RADseq to sequence thousands of markers in painted turtles from seven locations across the western range of the species. We found that populations were genetically differentiated, with patterns of divergence consistent with geographic distance among populations. Demographic analysis suggested that western populations of painted turtles are the result of serial founder effects during range expansion westward following the LGM, consistent with the decrease in genetic variation and increase in shorter ROH from eastern to western populations. This history of range expansion poses a serious challenge to accurate identification of loci under while the results from one outlier detection method can be explained by serial founder effects, another method found that many more outliers may have resulted from local adaptation and thus warrant further investigation. Together, we found candidate loci that are putatively under selection or are the result of genetic surfing during range expansion, and we conclude that serial founder effects have shaped genomic divergence among western populations of painted turtles.
In many widespread species, populations are characterized by strong genetic differentiation from the co-occurring forces of selection and genetic drift (e.g., Rödin-Mörch et al., 2019). During range expansions, and in the absence of significant gene flow, reduction in genetic diversity is often observed at range edges due to serial founder effects (DeGiorgio et al., 2011). Spatial bottlenecks, which may be expected in turtles that utilize freshwater habitat that is not contiguous across the western United States, can drive patterns of stronger genetic structure in populations at range edges than those from the original range (Excoffier et al., 2009; e.g., Graciá et al., 2013). Further, these bottlenecks increase the probability of genetic surfing, which can easily be misinterpreted as selection along clines in environmental conditions along the expansion front (Hoban et al., 2016). Thus, before assessing signatures of selection across the genome, we first investigated the population genomic structure of the western range of painted turtles to test past biogeographic hypotheses of the species’ range in North America.
Previous analyses of population genetic differentiation among western painted turtle populations have consistently revealed low genetic divergence. Starkey et al. (2003) investigated the phylogeography of painted turtles across the species’ range using a single mitochondrial marker (control region CR) and reported pairwise divergence of 0.151–1.368% within the western populations. Over 60 individuals sampled had the same haplotype, with 28 individuals sharing haplotypes with one or two differences from this main haplotype (Starkey et al., 2003). The authors concluded that this lack of divergence suggests a recent radiation of painted turtles into the western United States and Canada, despite fossil evidence suggesting their presence in the Great Plains as long ago as 1.9 million years (Holman, 1995). A more recent study expanded sampling to include a nuclear marker, PAX-P1, and found that the western range of painted turtles also included low haplotype diversity, though populations did not distinctly group together (Jensen et al., 2015). Finally, Reid et al. (2019) included eleven microsatellite loci to further elucidate population differentiation. Although fewer populations were sampled west of the Mississippi River (N=3), the results were similar in pattern to those presented here; populations were arranged from east to west in a principal coordinates analysis, and there was a significant relationship between FST and geographic distance (r = 0.85, Reid et al., 2019). Ecological niche modeling and demographic modeling suggested a single glacial refugium followed by population expansion (Reid et al., 2019), in contrast to Jensen et al. (2015), which did not find evidence of population expansion.
Despite the low genetic diversity found in the western range by previous studies, our study genetically distinguished most populations except turtles from Oregon and Idaho. Indeed, the FST values found among populations were high in some cases, particularly when comparing northwestern and New Mexico individuals (FST = 0.31). Given the filtering of SNPs in this study is conservative, and filters such as the HWE filter can remove loci that are strongly diverged, these FST values are likely underestimates (e.g., Pearman et al., 2022). These populations reside at the range edges for painted turtles, and as such support the pattern of stronger genetic divergence at range edges compared to the range center (Excoffier et al., 2009). Importantly, we detected patterns of genetic variation largely consistent with serial founder events driving population genomic structure. The evidence of timing and severity of bottlenecks in the western populations, with Idaho turtles having experienced the most recent and severe bottleneck following the LGM, and the decreasing genetic variation with distance from eastern populations both support range expansion from east to west. Further, the number of short ROH increased, consistent with bottlenecks leading to the current low levels of genetic variation in Oregon and Idaho, and to a lesser extent in the New Mexico population. ROH analysis with RADseq data is only an approximation of ROH dynamics across the genome, and the reference we used (from the northwestern US) may have slightly biased our ROH to be lower in eastern populations (by inflating heterozygosity, Thorburn et al., 2023). However, the number of markers used and the agreement with demographic results, which have been found to be less sensitive to reference genetic distance when divergence is less than 3% between reference and target species (Prasad et al., 2022), suggest that the increasing number of short ROH from eastern to western populations is likely a close representation of ROH in western painted turtles (Shafer et al., 2016). The conclusion of Reid et al. (2019) that the species’ widespread range is the result of the spread of populations westward after the retraction of glaciers following the late Pleistocene is confirmed comprehensively by the genome-wide sequencing performed here. In sum, our results suggest that serial founder effects during range expansion largely shaped genetic variation in western painted turtles.
Serial founder effects characterize not only the genetic structure of painted turtles, but also of other Testudines that experienced post-glacial colonization of northern ranges. For example, in the spur-thighed tortoise, populations in the northern range were characterized by reduced genetic diversity, strong population differentiation, and clinal variation from southern to northern populations in allele frequencies, all patterns consistent with genetic surfing (Graciá et al., 2013). Similar patterns were reported at microsatellite loci in the European pond turtle (Pereira et al., 2018).
We also find evidence consistent with multiple instances of human-mediated translocation of painted turtles from one population to another. We found two admixed Oregon individuals which appear to be the result of admixture between a northwestern turtle and one sourced from a Midwestern population. Individuals sequenced in a previous study also indicated potential evidence of human-mediated dispersal into the western range (Jensen et al., 2015). This result highlights the importance of managing and reducing the release of pet turtles into native populations, which remains a common practice and has led to both the establishment of invasive turtle populations and the spread of pathogens globally (e.g., Héritier et al., 2017). This concern was raised in the conservation assessment for the western painted turtle in Oregon (Gervais et al., 2009), and we can now establish that introductions have occurred even in areas that previously were not suspected of introductions in Oregon.
Another goal of this study was to evaluate genomic evidence of selection in western populations of painted turtles, which might be expected given the difference in climatic conditions experienced and the phenotypic variation among populations. Examples of phenotypic variation in these western locations include mean adult female body size and average clutch sizes, which vary latitudinally (Iverson & Smith, 1993), and the amount of patterning on the plastron and scute alignment patterns (Ultsch et al., 2001). Additionally, hatchlings from these sampled locations raised in common-garden conditions exhibited population-specific incubation durations and hatchling body mass (Bodensteiner et al., 2019). When we assessed genes that we a priori hypothesized would be under different selective pressures in the widespread sampling locations, we found little variation in genes related to TSD or anoxia tolerance across individuals that would support local adaptation for these traits. BayeScan found only 53 outlier loci (out of 212,670 SNPs), all of which had allelic patterns consistent with either selection or genetic surfing from east to west, as the SNPs were segregating in the eastern populations and fixed in New Mexico or Oregon and Idaho. This ambiguous finding is unsurprising, as BayeScan assumes that gene frequencies can be approximated by a Dirichlet distribution, and this assumption is violated under scenarios of range expansion (Foll & Gaggiotti, 2008; Lotterhos & Whitlock, 2014). Further, differences in allele frequencies from genetic surfing during range expansion can be indistinguishable from selection, and thus a history of range expansion can greatly hinder the ability to detect local adaptation in selection scans (Hoban et al., 2016; Hofer et al., 2009). Indeed, all FST outlier tests perform poorly, with many false positives, when range expansion shaped population genetic structure (Lotterhos & Whitlock, 2014; Luu et al., 2017). Our BayeScan outlier loci are similar to those found in invasive house finch populations, which could also be explained by genetic drift from founder events alone, rather than selection (Shultz et al., 2016).
Given the difficulties of detecting outlier loci in species with a history of range expansion, we used an additional method of outlier detection, pcadapt, which does not have the same assumptions with respect to demographic history as BayeScan. pcadapt instead detects outliers by finding loci that are highly correlated with ordination axes of the PCA rather than FST-related statistics (Privé et al., 2020). Analysis with pcadapt resulted in 228 outliers excluding HWE-filtered loci and 476 including them, many of which had allelic patterns that did not follow the predicted east-to-west pattern and yet were still strongly differentiated among populations. Additionally, one of the nonsynonymous mutations appears to be a de novo mutation in the New Mexico population that has risen to high frequency. The gene, ZAR1l, is important for coordination of mRNA translation during early embryonic development and is conserved across vertebrate lineages (Heim et al., 2022; Sangiorgio et al., 2008). Some of these outliers, particularly those that appear to have arisen de novo in specific populations, have reached high frequency, and were not segregating in the eastern sampling locations, may be the result of selection as painted turtles spread to inhabit the diverse environmental conditions of the western United States. However, we note that while pcadapt has lower FDR in range expansion scenarios than BayeScan, it is also prone to false positives in species with history of range expansion (Luu et al., 2017), and these loci should be studied further with alternative methods.
In species with long generation times, it can be difficult to assess local adaptation with reciprocal transplant experiments. Genomic methods can provide candidate loci for further investigation, but in cases of range expansion, it is clearly difficult to distinguish genetic surfing from selection. For example, in Alpine ibex, a species that experienced a recent bottleneck due to overharvesting, the influence of bottlenecks on genetic variation largely prevented the accurate identification of loci under selection (Leigh et al., 2021). Further, in geographically isolated populations with low levels of gene flow, we would anticipate that local adaptation of a trait would be achieved through divergence at very few large-effect loci and many more small-effect loci (Savolainen et al., 2013). Small-effect loci are often missed in FST outlier tests (Kemper et al., 2014), and thus we may not be able to effectively identify evidence of local adaptation in the phenotypes that vary among populations using these methods. Finally, in this study we used reduced-representation sequencing, which is an excellent tool for understanding population genetic structure (Catchen et al., 2017), but may miss regions of large differentiation consistent with phenotypic adaptation. Whole-genome sequencing may provide increased resolution for loci of large effect, and genome-wide association studies may be a better method for addressing local adaptation in complex phenotypes in species with a history of serial founder effects (Berg & Coop, 2014).
We studied the painted turtle, a widespread ectotherm, to further our understanding of how certain species inhabit entire continents. One frequent feature of widespread temperate ectotherms is a history of isolation in refugia, followed by postglacial colonization (Shafer et al., 2010). Indeed, we found serial founder effects consistent with range expansion following the LGM in this study. While local adaptation can facilitate colonization, long-lived, slowly evolving lineages may require much longer to locally adapt in order to facilitate postglacial spread. Alternatively, phenotypes that ensure survival in a wide range of environmental conditions may predate glaciation. In ectotherms, cold tolerance may be linked to achieving a widespread range during postglacial range expansion (Horreo & Fitze, 2018). Painted turtles possess multiple exceptional phenotypes to survive cold temperatures (e.g., anoxia tolerance and supercooling tolerance of hatchlings). These phenotypes are present across latitudes, suggesting they were acquired at least in some form before range expansion (Packard & Packard, 2003; Reese et al., 2004). Thus, long-lived, widespread temperate ectotherms may achieve their large ranges through pre-existing phenotypes. More generally, phenotypic plasticity is important for facilitating survival across environments during rapid colonization (reviewed in Schmid et al., 2019), particularly in ectotherms with TSD (e.g., Bodensteiner et al., 2023). Future work should continue to investigate the role of pre-existing phenotypes, local adaptation, and phenotypic plasticity in facilitating the large geographic ranges seen in some ectotherms. In conclusion, we clarified past demography and investigated signatures of selection in a long-lived, widespread ectotherm that faces an uncertain future as the global climate continues to rapidly change.