Discussion
The aim of the present study was to investigate genome-wide patterns of adaptive variation associated with climate in european bank vole populations. Searching for potentially adaptive loci in a multivariate framework is a powerful approach, especially because many adaptive traits are likely to be polygenic in nature (Barghi et al., 2020; Wellenreuther & Hansson, 2016)⁠. We used ddRAD sequencing and a combination of landscape genomic approaches to uncover environment-related evolutionary processes in the bank vole. We characterized adaptive genetic variation by using univariate GEA methods to detect outlier loci correlated with climate. We identified 74 genes of interest and functional annotation suggested that energy homeostasis and response to pathogen infection are important targets of spatially varying selection in the bank vole. In addition, we have shown that both population structure and climate play important (and common) roles in explaining neutral genetic differentiation across the bank vole range. Genetic variation among candidate loci was mainly associated with variation in annual mean temperature, highlighting the importance of this climatic variable in bank vole adaptation.
Loci with true selection signals must be distinguished from loci that exhibit genetic differentiation between populations caused by neutral forces. Correcting for these effects is an important concern when identifying candidate loci subject to selection. Proper correction can help to avoid possible spurious detection of candidate loci whose allele frequencies resemble signals of selection but are the result of neutral processes due to the shared history of populations (de Villemereuil et al., 2014). We therefore used GEA methods that correct for such confounding effects by accounting for population genetic structure in the bank vole.
Agreement between the different univariate methods was reasonable (14.9% and 15.1% of SNPs were detected by both methods for PC1 and PC2, respectively) and was consistent with results observed in other empirical studies using univariate methods (Harrison et al., 2017; Prates et al., 2018; Pritchard et al., 2018)⁠. The overlap between univariate methods and RDA was relatively small, with only 5.1% of the loci detected by RDA being also detected by LFMM or Bayenv2. A simulation-based study that tested the performance of univariate and multivariate GEA methods showed that the performance of these approaches varied depending on the strength of selection (Forester et al., 2018)⁠ and that RDA may be more robust to our sampling design that does not maximize environmental differentiation. The same study also suggests that combining results from univariate and multivariate approaches may help to increase power and to reduce false-positive rates. Our study supports these findings and also provides a strong argument for using multiple approaches when searching for signals of local adaptation in highly structured populations.
Population structure explained a large part of genomic variation, resulting in a strong pattern of isolation by distance. Even after taking into account the effects of climate variation, population structure still accounted for 33% of the total genomic variation explained by the RDA. These results are not surprising, as bank vole populations experience recurring population crashes and effective gene flow between populations is generally low, which results in isolation by distance across both smaller and larger geographic scales (Aars et al., 1998; Gerlach & Musolf, 2000; Guivier et al., 2011; Redeker et al., 2006), and although populations have strongly diverged in space since the last glaciation (for review Kotlik et al., 2022))⁠. On the other hand, climate explained 18.3% of total genomic variation when taking population structure into account, indicating association between bank vole genetic variation and environmental gradients. A large proportion of genetic variation (48.7%) could not be attributed either to the effects of climate or to spatial population structure alone, indicating that a large proportion of genomic variation associated with climate is geographically structured. The phylogeography of C. glareolus is marked by distinct genetic lineages, which resulted from survival within glacial refugia and recolonization of Europe at the end of the last glaciation (Filipi et al., 2015; Horniková et al., 2021; Kotlik et al., 2006; Marková et al., 2020). Post-glacial expansion from glacial refugia may result in clines of neutral allele frequencies coinciding with climate variables related to geography (Lotterhos & Whitlock, 2014; Rellstab et al., 2015)⁠, possibly explaining the large proportion of genomic variation. Our results suggest that annual mean temperature is an important driver of adaptive genomic variation and thus may be an important selection pressure influencing adaptation in the bank vole populations (Tiffin & Ross-Ibarra, 2014)⁠ as reflected by the strong association between temperature and polygenic scores. The latter implies that different alleles are maintained in different thermal environments, suggesting the presence of spatially varying selection pressure.
Temperature is one of the most important environmental factors affecting physiological processes such as the aerobic scope (Pörtner, 2001)⁠ and metabolism (Lovegrove, 2003)⁠, which in turn affect a variety of life history traits (Simons et al., 2011; Tökölyi et al., 2014)⁠. Numerous studies have associated clinal temperature variation and genome scans and found signals of selection in genes related to energy homeostasis and metabolism in endotherms (e.g., Andrew et al., 2018; Fumagalli et al., 2015; Hancock et al., 2011; S. E. Harris & Munshi-South, 2017; Harrison et al., 2017; Lv et al., 2014)⁠. This suggests that temperature is one of the most important environmental variables driving local adaptation. Indeed, temperature has been linked to adaptive genetic variation in other small mammals, such as populations of the recently introduced house mouse (M. musculus ) along a latitudinal cline in eastern North America (Phifer-Rixey et al., 2018)⁠ and populations of the climate-sensitive American pika (Ochotona princeps ) along an altitudinal cline (Waterhouse et al., 2018)⁠. The distribution and abundance of C. glareolus from the Eastern lineage in a contact zone with the Carpathian lineages correlated negatively with July temperature, suggesting that these individuals are better adapted to cooler conditions (Tarnowska et al., 2016)⁠, supporting temperature as a driver of adaptive genetic variation in the bank vole.
Artificial selection experiments for higher aerobic exercise performance in bank voles resulted in an increase in resting metabolic rate and thus resulted in the development of increased cold tolerance (Sadowska et al., 2015; Stawski et al., 2017)⁠. This argues for a genetic basis for thermal adaptation in bank voles that may allow individuals under natural conditions to adapt to colder environment by having more energy available for thermogenesis (Stawski et al., 2017)⁠. Although the selection regime increased cold tolerance, it also decreased the ability to thermoregulate at higher temperatures (Grosiak et al., 2020). This suggests that warmer temperatures may also be difficult for small mammals to cope with, as this can easily lead to overheating (Rezende et al., 2004). This in turn could also lead to specific metabolic adaptations in populations at warmer climates due to increased selection pressure. In this study, AMT differed between -1.4°C (NE3.fi) and 12.4°C (S.it) among populations. Thus, bank vole populations in Europe are exposed to different environmental temperatures, likely resulting in different energetic requirements and adaptive genetic divergence in metabolic traits throughout the species’ range.
We have identified a number of promising candidate genes that could be considered for future research aimed at linking phenotypic and genotypic variation. The function of these candidate genes provides insight into the physiological processes that may have experienced selection across climatic gradients. Different populations are exposed to different local climatic conditions that determine the energy requirements, diet and different pathogen or predator communities. In this context, we have identified a number of candidate genes related to lipid metabolism and the immune system that appear to be subject to temperature-related selection.
Adipose tissue plays an important role in energy homeostasis and accounts for a large portion of the energy reserves of small mammals (Birsoy et al., 2013; Sethi & Vidal-Puig, 2007)⁠. In particular, brown adipose tissue is important for metabolic heat production through non-shivering thermogenesis under cold conditions (Cannon & Nedergaard, 2004; Klaus et al., 1988). Two candidate genes associated with lipid metabolism are therefore of particular interest : the PRIP gene encodes an enzyme that modulates lipid metabolism and serves as a signalling molecule for non-shivering thermogenesis (Kanematsu et al., 2019; Oue et al., 2016), and LRRC8C which encodes a structural component of the volume-regulated anion channel in adipocytes and is associated with the early phase of adipocyte differentiation and diet-induced obesity (Hayashi et al., 2011; Tominaga et al., 2004)⁠.
Other candidate genes with functions related to energy homeostasis include NTRK2 , which encodes the TrkB-receptor critical for maintaining energy homeostasis by controlling food intake and body weight and is responsible for regulating adaptive thermogenesis (Houtz et al., 2021; Xu & Xie, 2016). Finally, the product of IGF1 has wide ranging effects on metabolism by coordinating protein, carbohydrate, and lipid metabolism in a variety of different cell types (Baker et al., 1993; Laron, 2001)⁠. Several of the candidate genes are associated with obesity in humans including DNAH8 (Söhle et al., 2012)⁠, IGF1 (Berryman et al., 2013)⁠, KCNH1 (Vasconcelos et al., 2016)⁠, LRRC8C (Hayashi et al., 2011), NTRK2 (Gray et al., 2007)⁠ and PRIP (Yamawaki et al., 2017)⁠, suggesting that they play a role in controlling energy homeostasis.
The results showed significant enrichment of genes related to the regulation of respiratory burst. The respiratory burst plays an important role in the immune system. It is a crucial reaction that occurs in phagocytes to degrade internalised pathogens after phagocytosis (Iles & Forman, 2002). In this context, we have identified a number of candidate genes that play important roles in the immune system. For example, the product of DUSP10 , which was associated with this significant GO term, plays an important role in regulating both innate and adaptive immune responses through its regulatory influence on the MAPK pathway (Arthur & Ley, 2013; Seternes et al., 2019). Two other candidate genes, BATF3 and BACH2 , both encode transcription factors that regulate T helper cell function. Interestingly, they also interact with each other to bind to regulatory regions of cytokine gene loci and prevent excessive T helper response (Kuwahara et al., 2016; Yamashita & Kuwahara, 2018)⁠. Another candidate gene of interest is STAT4 . This gene encodes a transcription factor responsible for the differentiation of T helper cells (Kaplan, 2005)and is part of the JAK-STAT signalling pathway that controls the immune response to viral infections (Villarino et al., 2017)⁠. JAK-STAT is one of the significantly enriched signaling pathways associated with Puumala hantavirus infection in the bank vole (Rohfritsch et al., 2018).  It has also been found to play a role in the immune response to Sin Nombre hantavirus in deer mice (Peromyscus maniculatus ) (Schountz et al., 2012, 2014)⁠. This suggests that alterations in this gene may be related to Puumala hantavirus infections in the bank vole populations we studied. Similar evidence for differential selection on immune-related genes has been observed in bank vole populations along environmental gradients at both broad and local scales, using candidate genes (Dubois et al., 2017; Guivier et al., 2014) and more exploratory genome-wide approaches (Rohfritsch et al., 2018; White et al., 2013)⁠. For humans, the diversity of the local pathogenic environment is the predominant driver of local adaptation (Fumagalli et al., 2011) and we may speculate that pathogen environment and pathogen pressure is closely associated to climate variation, with rapid adaptations expected due to climate change. Taken together, the selection signals in the most promising candidate genes (Table 2, full overview in Table S5) suggest that the energy balance and immune system in the bank vole are important targets of temperature-mediated, spatially varying selection.