Adaptive associations with the environment
As we hypothesized, some populations within each species exhibited signatures of adaptation to local environmental conditions. At the interspecific level, SIDGS were generally more associated with grassland and higher soil particle size, which is positively correlated with soil temperature, while associations of NIDGS were varied and appeared to vary mostly on a population basis (Figures 2 and S3A, Supporting information). These differences are within our expectations of adaptive differentiation between the two species, given that NIDGS and SIDGS occur in high and low elevations, respectively (Hoisington-Lopezet al. 2012). These genomic signatures of adaptive differences between NIDGS and SIDGS might also reflect associations with timing of snow melt, or site productivity, which are related to difference in hibernation (Hoisington-Lopez et al. 2012; Zero et al.2017; Goldberg, Conway, Mack, et al. 2020).
Intraspecific differences can quickly arise in populations that become isolated and/or that are distributed across highly variable environmental gradients (Doebeli & Dieckmann 2003; Smith et al.2019), and there is an important role of adaptive variants and their specific adaptations in the maintenance of genetic diversity and the long-term persistence of threatened species (Rubidge et al.2012). Small populations can naturally become locally adapted in the process of isolation, resulting in a gradual increase in adaptive differentiation from other populations, compared to neutral differentiation (Doebeli & Dieckmann 2003; Holderegger et al.2006; Wood et al. 2016). In agreement with our initial hypothesis (c), NIDGS showed a higher number of distinct populations both the neutral and adaptive level, with Rocky Top and Lower Butter as the most distinct populations and Tamarack to a lesser extent (Figures S11A and S11B, Supporting information). None showed particularly low genetic diversity or significantly lower H O thanH E (indicative of inbreeding), indicating that these populations are unique but also demographically stable, despite local bottlenecks (Assis et al. 2013). Rocky Top and Lower Butter were mainly associated with environmental variables relating to elevation (Figure 3), and they are the highest elevation populations sampled in this study (around 1700 and 1600 m above sea level, respectively). This suggests that these populations might have particular adaptations to elevation compared to other NIDGS populations, a pattern also observed in other North American ground squirrels (Eastman et al. 2012). Tamarack is located on a west facing slope which remains colder until later during the day, similarly to populations at higher elevation sites. Tamarack was also found to be associated with increased soil particle size (Figure 3), which is positively correlated with soil temperature (Figure S3B, Supporting information). Combined, Tamarack NIDGS may reflect an adaptation to the colder temperatures despite living at ~1,275m (relatively lower elevation for NIDGS total range). The highest adaptive differentiation in SIDGS was observed for Paddock (PA), with patterns of local adaptation particularly associated with grassland land cover rather than shrub/scrubland, as well as increase Annual Mean Temperature and Isothermality. Greater densities of SIDGS have been associated with increased cover of perennial grasses and diversity of native perennial plants (Lohr et al. 2013). Additionally, burrow density appears to be influenced by the presence of native forbs typical of shrub/scrubland, which was an environmental variable associated with adaptive variation of the remaining SIDGS populations (Figure 6C). Forbs have a large influence in protein intake for fat storage before hibernation, and thus it is a strong predictor of overwintering survival of all age classes, even when present at low densities (Barrett 2005). The results of the GEA could indicate a particular adaptation of Paddock to areas with lower forb density, but further sampling of other SIDGS populations, as well as diet studies would be needed to confirm this hypothesis.
We did not find particular Gene Ontology terms associated with putative functions of genes in our analysis, which could be related to a lack of annotations of the thirteen-lined ground squirrel genome, or also the fact that no particular Gene Ontology term stands out as causative of the differentiation of the two species (Szkiba et al. 2014). Still, we found two non-synonymous substitutions in known genes identified by the RDA, one of which was also identified bypcadapt , the NPR1. This gene has been found to be highly expressed in brown adipose tissue of thirteen-lined ground squirrels during hibernation, relating to heat production during periodic arousals from hypothermic torpor (Hampton et al. 2013). Adaptations relating to adipose tissue are common in high altitude adapted rodents (Gossmann et al. 2019), and these associations suggest that there are adaptive differences between NIDGS and SIDGS relating to their hibernation patterns. In fact, one study on Columbian ground squirrel (Urocitellus columbianus ) found that emergence timing was heritable and hence, squirrels’ emergence date may also reflect local adaptations (Lane et al. 2011). In accordance to our hypothesis (d), altitude differences might have resulted in adaptive differences between NIDGS and SIDGS in production and storage of fat, and increased metabolism and oxygen transport at higher elevations, as seen in other mammals (Faherty et al. 2018; Waterhouse et al. 2018; Garcia-Elfring et al. 2019).