Gene regulation underlies environmental adaptation in house mice.

Changes in cis-regulatory regions are thought to play a major role in the genetic basis of adaptation. However, few studies have linked cis-regulatory variation with adaptation in natural populations. Here, using a combination of exome and RNA-seq data, we performed expression quantitative trait locus (eQTL) mapping and allele-specific expression analyses to study the genetic architecture of regulatory variation in wild house mice (Mus musculus domesticus) using individuals from five populations collected along a latitudinal cline in eastern North America. Mice in this transect showed clinal patterns of variation in several traits, including body mass. Mice were larger in more northern latitudes, in accordance with Bergmann's rule. We identified 17 genes where cis-eQTLs were clinal outliers and for which expression level was correlated with latitude. Among these clinal outliers, we identified two genes (Adam17 and Bcat2) with cis-eQTLs that were associated with adaptive body mass variation and for which expression is correlated with body mass both within and between populations. Finally, we performed a weighted gene co-expression network analysis (WGCNA) to identify expression modules associated with measures of body size variation in these mice. These findings demonstrate the power of combining gene expression data with scans for selection to identify genes involved in adaptive phenotypic evolution, and also provide strong evidence for cis-regulatory elements as essential loci of environmental adaptation in natural populations.

Understanding the genetic basis of adaptation is a major goal in evolutionary biology. Cis-regulatory mutations, which can change the expression of proximal genes, have long been predicted to be important targets for adaptive phenotypic evolution (King and Wilson 1975; Wray 2007; Stern and Orgogozo 2008; Wittkopp and Kalay 2012). One reason for this is that cis-regulatory mutations may have fewer deleterious pleiotropic effects than protein-coding changes. While protein-coding mutations may affect protein products across tissues and developmental stages, cis-regulatory mutations can affect the expression of genes in spatially and temporally specific ways. In apparent support of this idea, several studies have identified positive selection on noncoding regions (e.g., Jenkins et al. 1995; Crawford et al. 1999; Kohn et al. 2004; Andolfatto 2005; MacDonald and Long 2005; Holloway et al. 2007; Jeong et al. 2008; Torgerson et al. 2009) and an important role for noncoding variation in local adaptation (e.g., Jones et al. 2012; Fraser 2013).

Despite the accumulating evidence that regulatory loci play an important role in adaptive evolution, there are still only a handful of cases where cis-regulatory mutations have been linked to ecologically important traits. Among the best examples are adaptive coat color differences in deer mice (Linnen et al. 2013), the ability to digest lactose in humans (Tishkoff et al. 2007), and pelvic reduction in sticklebacks (Chan et al. 2010). Most examples of adaptive gene expression have been identified through candidate gene approaches, which typically favor traits for which components of a pathway are already known and the genetic basis of the trait is relatively simple. However, most traits are influenced by many loci of small to modest effect. Thus, identifying genetic variants associated with adaptation at complex traits is key to understanding the genetic basis of adaptation.

One avenue for linking adaptive non-coding variation to either molecular or organismal phenotypes is through gene expression. In expression quantitative trait loci (eQTLs) mapping, gene expression levels are tested for associations with genetic markers to identify variants that contribute to expression phenotypes. Expression quantitative trait mapping is an effective method for identifying regulatory variants because gene expression is frequently influenced by nearby cis-eQTLs (Nica and Dermitzakis 2013). Cis-eQTLs have been successfully detected with small sample sizes (Montgomery and Dermitzakis 2011; Tung et al. 2015) and in wild individuals from natural populations (Tung et al. 2015). Combining eQTL mapping with genomic scans for selection can be a powerful method for identifying the gene targets of adaptive genetic variation (Fraser 2013; Ye et al. 2013) and potentially linking this variation to adaptive organismal phenotypes.

House mice (Mus musculus domesticus) provide a useful model for studying the genetic basis of adaptation. House mice are an important biomedical model and have a distribution that mirrors that of human populations (Phifer-Rixey and Nachman 2015). In the eastern United States, house mice show latitudinal variation consistent with local adaptation. Mice collected at northern latitudes are heavier than mice at southern latitudes, and their progeny also show differences in a common laboratory environment, indicating that this difference is genetic (Lynch 1992; Phifer-Rixey et al. 2018). This observation conforms to the classic ecogeographic observation known as Bergmann's rule that animals in colder climates have larger mass to reduce heat loss (Bergmann 1847). While Bergmann's rule has been observed in many groups, including humans (Ashton et al. 2000; Ruff 2002; Foster and Collard 2013), no study so far has linked this pattern to variation at specific genes. Consistent with energetic adaptation of mice from eastern North America, laboratory strains founded from northern and southern locations also show differences in aspects of blood chemistry, including leptin, glucose, and triglyceride levels (Phifer-Rixey et al. 2018).

Recent work with these populations identified hundreds of genes with environmental associations in North American (Phifer-Rixey et al. 2018). Here, we combine a genomic scan for selection with expression quantitative trait mapping to identify regulatory variants that contribute to gene expression differences and show signals of selection in these populations, identifying two strong candidate genes for adaptive phenotypic variation. To our knowledge, this study represents the first case where genomic scans have been combined with eQTL mapping to identify regulatory variants in natural populations that underlie an adaptive organismal phenotype.

Results

Cis-regulatory variation in wild house mice

To characterize regulatory variation in wild mice, we sequenced liver transcriptomes from 50 mice collected from five populations along a latitudinal transect on the east coast of North America (Fig. 1; Supplemental Table S1; Supplemental File S1). Mice were collected from 29°N to 44°N latitude. Liver was collected in RNAlater, and body mass and length were recorded for each individual. From these individuals, we produced a total of ∼1.2 billion RNA-seq reads with an average of 15,473,949 uniquely mapped exonic reads per sample, which were used to quantify gene-wise mRNA abundance (hereafter, gene expression). We also analyzed DNA sequence data generated from exome-capture of the same individuals (Phifer-Rixey et al. 2018). Exome and RNA-seq data were used to identify variants segregating in M. m. domesticus (see Methods).

Figure 1.

(A) Sampling locations along the east coast of North America (climate map obtained from NOAA, National Weather Service). (B) Consistent with Bergmann's rule, body mass in mice increases with increasing latitude (Pearson's correlation = 0.34, P = 0.018) (see Supplemental Table S9; Supplemental Methods).

We identified cis-regulatory variation using two complementary approaches, expression quantitative trait loci mapping and allele-specific expression (ASE). To identify cis-eQTLs, we tested for associations between variants within 200 kb of a gene and expression level using a linear mixed model. Variants near a gene are more likely to act in cis to affect gene expression. Cis-eQTLs typically have larger effect sizes than trans-eQTLs, making them easier to detect in small sample sizes (Montgomery and Dermitzakis 2011). After filtering, a total of 406,999 variants were identified using exome data and tested for associations with expression at 13,080 genes. We identified cis-eQTLs for 849 of these genes (6.5% of genes surveyed). Reflecting the probe set, the majority of cis-eQTLs were identified in gene bodies (57%) and introns (18%) (Supplemental Fig. S1).

Allele-specific expression (i.e., differences in expression between parental alleles) can also be used to infer epigenetic or genetic variation acting in cis (Cowles et al. 2002). As the two parental alleles are exposed to the same trans-acting environment within an individual, differences in expression at heterozygous sites can be used to infer cis-regulatory variation. A total of 28,234 exonic heterozygous sites, corresponding to 6738 genes, could be tested for ASE. Across all individuals, we found evidence for ASE for 442 genes at a false discovery rate of 5% (6.7% of genes surveyed) (Supplemental Table S2).

In investigating the power to detect cis-regulatory variation, we found that cis-eQTLs were more likely to be detected when SNP density is higher near and within the gene of interest (Mann-Whitney U test, P < 2.2 × 10−16) (Supplemental Fig. S2). We were more likely to detect ASE for genes with higher expression and higher SNP density (Mann-Whitney U test, P = 3.1 × 10−11 and P < 2.2 × 10−16, respectively) (Supplemental Figs. S2, S3). While differences in the power to detect ASE and cis-eQTLs can lead to the identification of different gene sets, we found significant overlap between the gene sets identified with these analyses (hypergeometric test, P = 5 × 10−6) (Supplemental Table S2).

Evidence for adaptive regulatory variation

To assess whether the regulatory variation documented above underlies adaptive differences among populations, we studied sequence and gene expression variation along a latitudinal cline (Fig. 1A). Clinal patterns of variation can reflect local adaptation as a response to spatially varying selection (Endler 1977). Regulatory variants with clinal frequencies that mediate clinal patterns of gene expression would be strong candidates for adaptive regulatory evolution. To identify such variants, we searched for cases where (1) gene expression is clinal, (2) gene expression is associated with a cis-eQTL, and (3) allele frequencies of the cis-eQTL vary clinally (Fig. 2). While geographic clines may alternatively be explained by isolation by distance, there is no evidence for isolation by distance for these populations (see Supplemental Methods).

Figure 2.

Overlap between genomic scans identifies regulatory variants that are candidates for clinal adaptation. (A) The LFMM |z-scores| for each SNP vs. chromosome position. SNPs with |z-scores| > 2 were considered clinal outliers. (B) Manhattan plot of cis-eQTL. Shown in red are significant SNPs. (C) Manhattan plot of gene starting position versus the correlation between gene expression and latitude. Points labeled in orange are genes for which expression is significantly correlated with latitude (P < 0.05). On the outside are ideograms with the location of genes for which these three signals (A–C) overlap. Figure created with Circos (Krzywinski et al. 2009).

To identify clinal patterns of gene expression, we tested for correlations between latitude and expression levels in the liver transcriptomes of the 50 wild individuals. We identified 1488 genes for which expression was significantly correlated with latitude (P < 0.05), 132 of which were associated with a cis-eQTL (Fig. 2). We also tested for differential expression between the most northern population (New Hampshire/Vermont) and the most southern population (Florida) and identified 458 genes with differential expression between the ends of the cline (Supplemental Fig. S4), 48 of which were associated with a cis-eQTL (at q < 0.1) (Supplemental Table S3).

To connect these patterns to clinal sequence variation, a genome scan using the program Latent Factor Mixed Models (LFMM) was performed to test for correlations between latitude and genetic variation while accounting for population structure (see Methods; Frichot et al. 2013). For this study, LFMM has an advantage over other methods because it does not assume a specific demographic model but still accounts for demographic history by estimating genome-wide covariance among allele frequencies. We focused on SNPs in the 5% tail of the distribution and considered these clinal outliers (|z-scores| > 2) (Fig. 2A). Blocks of linkage disequilibrium (LD) (Supplemental Fig. S5; Gabriel et al. 2002) were then inferred to identify colocalization between outlier SNPs and cis-eQTLs. Of cis-eQTLs that fell within the same LD block as an outlier, 17 were associated with genes that also show significant clinal patterns of gene expression (Table 1; Fig. 2; Supplemental Table S4). When comparing the latitudinal extremes, average estimates of F for these candidate loci were significantly higher than that of the full list of loci (full list average F = 0.10, candidate average F = 0.34; permutation test, P = 0.0014). Eight of these genes were also significantly differentially expressed between the ends of the cline (Supplemental Table S5). These 17 genes represent cases where cis-eQTLs contribute to expression differences between populations and show signals of local adaptation, making them strong candidates for adaptive regulatory variation.

Table 1.

cis-eQTLs that colocalize or are within the same LD block as a clinal outlier that also show expression changes correlated with latitude

Linking adaptive regulatory variation to specific traits

The liver plays a central role in metabolic processes in the body, and regulatory changes in this tissue may contribute to latitudinal variation in traits related to metabolism. Body mass varies clinally (Fig. 1A), and lab-born progeny from populations at the ends of the transect also show differences in blood glucose, triglyceride, adiponectin, and leptin levels (Lynch 1992; Phifer-Rixey et al. 2018). Four of the 17 candidate genes identified as strong candidates also have mutant phenotypes related to body weight and metabolism. Laboratory mutants for Cox7c and Hmgb1 are associated with changes in glucose levels (Blake et al. 2017), and mutants for Adam17 and Bcat2 are also associated with changes in body mass (Wu et al. 2004; She et al. 2007; Gelling et al. 2008; Blake et al. 2017), glucose (Serino et al. 2007; She et al. 2007; Blake et al. 2017), leptin (She et al. 2007; Gelling et al. 2008), and adiponectin levels (Serino et al. 2007; Blake et al. 2017). Another gene identified in this analysis, Iah1, transcriptionally regulates genes with important roles in lipid metabolism and triglyceride synthesis and falls under a QTL for fatty liver in mice (Kobayashi et al. 2016; Suzuki et al. 2016).

Adam17 and Bcat2 are candidates for adaptive differences in body mass

While knockout models can provide a link between genotypes and putative phenotypes, these models may not reflect the phenotypic consequences of mutations found in natural populations (Palopoli and Patel 1996). Changes in body weight are also among the most common effects of gene knockouts in mice and may often reflect downstream consequences of other phenotypic changes (Reed et al. 2008; White et al. 2013). While identifying the genetic basis of complex adaptive traits is challenging, gene expression provides an intermediate phenotype that may link sequence variants to organismal traits. To connect adaptive variation in body mass in these populations to genetic variation, we asked whether body mass differences were associated with gene expression differences in the set of candidate genes (Table 1). Since latitude and body mass covary in this sample (Fig. 1B), we controlled for latitude by regressing it out as a variable. We identified two genes, a disintegrin and metallopeptidase domain 17 (Adam17) (Fig. 3A–F) and branched chain amino acid transaminase 2 (Bcat2), for which expression was significantly correlated with body mass, after accounting for latitude as a covariable (Pearson's correlation, Adam17: P = 4.6 × 10−4, R2 = 0.22; Bcat2: P = 4.5 × 10−3, R2 = 0.17) (see also Supplemental Table S6; Supplemental Fig. S6). To further account for the possible confounding effects of population structure, we also looked at the correlation between expression level and body mass within each of the five populations. Replicating the pattern seen across populations, Adam17 expression was negatively associated with body mass in four of the five populations, and Bcat2 expression was positively associated with body mass in four of the five populations (Supplemental Figs. S7, S8). Despite a lack of power for within-population comparisons, the association between Adam17 expression and body mass was significant in New Hampshire/Vermont (Pearson's correlation, P = 3.5 × 10−3) and the association between Bcat2 expression and body mass was significant in Pennsylvania (Pearson's correlation, P = 0.03) and Georgia (Pearson's correlation, P = 1.8 × 10−3).

Figure 3.

Adam17 is a candidate for adaptive differences in body mass among mice in eastern North America. (A) Expression of Adam17 is correlated with latitude (P = 0.032, Pearson's correlation = −0.30). Sex was not a significant predictor of Adam17 expression. (B) A SNP at Chr 12: 21,332,631 was identified as a cis-eQTL for Adam17. (C) Allele frequencies of Chr 12: 21,332,631 in five populations. (D) The LFMM |z-scores| for sites on Chromosome 12 versus position. Points above the red line were considered clinal outliers in this study. The red box represents the peak in which Chr 12: 21,332,631 is found. (E) Nearby outlier SNPs in LD with Chr 12: 21,332,631. Correlations (r2, %) are given in each block. The z-scores for each site's association with latitude are given in parentheses. (F) Adam17 expression is significantly associated with body mass when controlling for latitude (Pearson's correlation, P = 4.6 × 10−4, R2 = 0.22). (G) Genotype at Chr 12: 21,332,631, the cis-eQTL for Adam17, significantly trends with body size when latitude is controlled for (Cochran-Armitage trend test, P = 0.034).

The cis-eQTLs for Adam17 and Bcat2 explain 34% and 29.7% of the variance in expression for these genes, respectively. Genotypes at these sites were also associated with differences in body mass (Mann-Whitney U test, Bcat2: TT > CC, P = 0.024; Adam17: CC > TT, P = 0.036) (Supplemental Fig. S9). Again, covariation between latitude and body mass can confound relationships between body mass and candidate genes. After regressing latitude from body mass to control for covariation between these variables, the Adam17 cis-eQTL was significantly associated with body mass (Cochran-Armitage trend test, P = 0.034) (Fig. 3G), although the Bcat2 cis-eQTL was not (Cochran-Armitage trend test, P = 0.14). The Adam17 and Bcat2 cis-eQTLs explain an estimated 8.35% and 1.51% of the variation in body mass, respectively. These estimates should be treated as approximations since they may be influenced by (1) unmeasured environmental differences between populations, (2) population structure (even when population structure is accounted for using principle components, as was done here) (see Browning and Browning 2011; Dandine-Roulland et al. 2016), (3) imperfect linkage disequilibrium between the surveyed SNPs and causal variants (Wray et al. 2013), and small sample size (Xu 2003). Nonetheless, it is likely that the effect size for the Adam17 cis-eQTL is large compared to what is seen in most human GWAS for complex traits (Stranger et al. 2011). Large-effect mutations may be favored in situations where populations are initially far from an optimum (Orr 1998; Dittmar et al. 2016). For example, variation at one gene accounts for a >2-kg weight difference between Europeans and Inuits (Fumagalli et al. 2015), and a single IGF1 allele in dogs accounts for 15% of variance in dog skeletal size (Sutter et al. 2007). House mice in this transect descended from mice in western Europe adapted to a Mediterranean climate and thus likely experienced strong selection pressures in a novel environment, potentially favoring some mutations of large effect.

To investigate regulatory variation at Adam17 in western Europe, we retrieved available liver RNA-seq and genomic data from European mice (Harr et al. 2016). We found that the Adam17 cis-eQTL is segregating within European populations (Supplemental Fig. S10A) and is significantly associated with liver expression in European individuals (P = 3.2 × 10−6) (Supplemental Fig. S10B; see Supplemental Methods). This suggests that adaptation by the large-effect regulatory variation at Adam17 in the United States is a product of selection on standing genetic variation.

Notably, Adam17 and Bcat2 are the two candidate genes from Table 1 with known lab mouse mutants that affect body mass (Wu et al. 2004; She et al. 2007; Gelling et al. 2008; Blake et al. 2017). Bcat2 encodes a protein that catalyzes the first step of branched-chain amino acid (BCAA) metabolism, which affects metabolism and body mass in humans and rodents (Newgard et al. 2009). Adam17 encodes a protein that regulates several signaling pathways. Adult Adam17 heterozygous and null mutants show differences in metabolic phenotypes including body mass, susceptibility to diet-induced obesity, and energy homeostasis (Serino et al. 2007; Gelling et al. 2008). ADAM17 and its physiological inhibitor, TIMP3, have also been reported to be involved in glucose homeostasis and adipose, hepatic, and vascular inflammation in both genetic and nutritional models of obesity in mice (Fiorentino et al. 2010; Menghini et al. 2012; Matsui et al. 2014). In addition to its association with body mass and metabolism in mice, in humans variation at ADAM17 has been linked to differences in body weight, BMI, waist circumference, and obesity risk (Junyent et al. 2010) and shows signatures of selection (Pickrell et al. 2009; Parnell et al. 2010; Fumagalli et al. 2011).

One target of ADAM17 activity is the epidermal growth factor receptor (EGFR) signaling pathway (Lee et al. 2003). Phenotypes observed in mice with mutant EGF receptors (including changes in body weight [Blake et al. 2017]) suggest that changes in EGFR signaling as a consequence of deficient ADAM17 activity may contribute to the metabolic phenotypes seen in Adam17 mutants (Gelling et al. 2008). We tested for an overrepresentation of genes in the EGFR signaling pathway in the set of genes with clinal expression by annotating genes to pathways using the PANTHER database (Thomas et al. 2003). We saw a 1.57-fold enrichment of genes in this pathway compared to a background set of genes expressed in the liver (hypergeometric test, P = 0.018). We also find that the gene that encodes the only known physiological inhibitor of Adam17, Timp3 (Le Gall et al. 2010), is differentially expressed between the northern and southern populations (q = 0.09) (Supplemental Fig. S11) and has expression that is correlated with that of Adam17 (Pearson's correlation, P = 0.02, R2 = 0.09) (Supplemental Fig. S11). Unlike Adam17, Timp3 expression is not associated with body mass (Pearson's correlation, P = 0.054), although our sample size may not be sufficient to detect an association.

The data above clearly suggest that regulatory variation at Adam17 and Bcat2 underlies adaptive differences in body mass, but they do not identify the specific causal mutations. To identify candidate casual mutations, we used annotations from the mouse ENCODE project (Mouse ENCODE Consortium et al. 2012) to search for putative regulatory elements near the Adam17 and Bcat2 cis-eQTLs. The Adam17 cis-eQTL is in LD with SNPs through a proximal enhancer and in the Adam17 promoter, both of which are active in the livers of adult mice. Low-coverage whole-genome data show that there are variants segregating within this enhancer in these populations (Supplemental Figs. S12, S13; whole genome data from Phifer-Rixey et al. 2018). Two of the Adam17 promoter variants are also clinal outliers (Supplemental Fig. S14). The Bcat2 cis-eQTL is within an intronic region and is not in LD with annotated regulatory elements that are active in liver tissue.

Expression modules are correlated with body size variation in natural populations of house mice

Finally, we used a gene co-expression network approach to identify biologically related gene sets associated with phenotypic variation in these populations. Weighted Gene Co-expression Network Analysis (WGCNA) was used to identify groups of genes with highly correlated expression, called co-expression modules (see Methods; Langfelder and Horvath 2008). Expression modules were assigned for male and female mice separately, and then male-female consensus modules were created to identify co-expression patterns shared across sexes.

Co-expression modules were then tested for correlations with measures of body size (Supplemental Figs. S15–S17). Five expression modules in males and five expression modules in females were correlated with trait variation (Supplemental Fig. S18). Trait-associated modules were enriched for a number of Gene Ontology (GO) categories compared to the background set of genes expressed in the liver, including growth factor binding (q = 5.3 × 10−8) and lipid metabolic process (q = 1.2 × 10−2). None of the male-female consensus modules were significantly correlated with organismal traits, indicating that associations between co-expression modules and traits are sex-specific (Supplemental Fig. S17).

Focusing on the modules with the highest trait correlations (royalblue module in females, corr = 0.92, P = 2 × 10−8 and black module in males, corr = 0.8, P = 5 × 10−8, for body mass index), we annotated genes with mutant phenotypes collected from Mouse Genome Informatics (MGI) (Blake et al. 2017). Supporting the association between these expression modules and phenotypic variation, we found that many of the genes with high connectivity in these modules have mutant phenotypes related to body size or metabolism (Fig. 4). For example, the most connected gene in the female royalblue module is Nr2c2. Mutant phenotypes for Nr2c2 include changes in eating behavior, energy homeostasis, body mass, size, and blood chemistry. Similarly, highly connected genes in the male black module (e.g., Col3a1, Col1a1, Col1a2, Col5a2, Sparc, Bcam, Fstl1, Igfbp5, Cpe, Cav1, Lamc1, Ltbp3, Krt7) show mutant phenotypes related to body mass and body size. Four of these genes (Adamts2, Col1a1, Col1a2, Sparc) were also identified as hub genes in the module most highly correlated with mouse body weight in another study utilizing an F2 laboratory cross (Ghazalpour et al. 2006).

Figure 4.

Visualization of the most connected genes in the female “royalblue” (A) and the male “black” (B) co-expression modules with VisANT (Hu et al. 2008). The royalblue module is associated with BMI (P = 2 × 10−8) and body length variation (P = 6 × 10−6). The black module is associated with BMI (P = 5 × 10−8), body mass (P = 0.001), and body length variation (P = 3 × 10−10). Blue circles represent genes for which we identified a cis-eQTL that explains a component of expression variation. Circles with black borders are genes with mutant phenotypes related to body size or metabolism. Phenotype information was collected from MGI (Blake et al. 2017).

Finally, we used the co-expression data set to identify regulatory variation within modules associated with body size. Within the body size-associated modules (Supplemental Fig. S18), we associated 189 genes with a cis-eQTL, including several highly connected genes in the sex-specific modules with the highest trait correlations (Fig. 4). As in the previous analysis, we then searched for genes with a cis-eQTL that colocalized with a clinal sequence variant. We identified 15 genes with clinally varying cis-eQTLs in the body size-associated modules (Supplemental Table S7). We found that gene expression for four of these 15 genes was significantly correlated with BMI in one sex (females: Ube2q2, P = 0.0002; 3110082I17Rik, P = 0.0027; Cep85, P = 0.017; males: Pygb P = 0.035). Cis-eQTLs associated with these genes were not significantly associated with BMI; however, our study is also underpowered for identifying sex-specific associations. The correlation between gene expression and BMI and the presence of clinal cis-eQTLs make these genes of interest for future study.

Discussion

Identifying loci and genes that underlie adaptive variation within and between populations is a major goal in evolutionary biology. One method used to identify such variants is genomic scans for selection. While many genomic scans attempt to link sequence variants to phenotypes through gene annotations and knockout models, most fail to connect genotypes to phenotypes in natural populations. Here, we used expression data from natural populations of house mice collected along an environmental gradient to link regulatory variation at two genes (Adam17 and Bcat2) with body mass variation. We have linked these genes to body mass variation by (1) associating cis-eQTLs with the expression of Adam17 and Bcat2, (2) associating the Adam17 and Bcat2 cis-eQTLs with body mass variation, and (3) associating the expression of these two genes with body mass variation. Supporting the association we see between these genes and body mass, mutant alleles for Adam17 and Bcat2 in laboratory mice are associated with changes in body mass and metabolism (Wu et al. 2004; Serino et al. 2007; She et al. 2007; Gelling et al. 2008; Blake et al. 2017). These two genes account for a substantial proportion of phenotypic variation in body mass among the mice studied here, with large effect sizes compared to those measured in GWAS for most complex traits. For traits under stabilizing selection within populations (as in virtually all human GWAS), effect sizes are expected to be much smaller than in comparisons between populations experiencing strong divergent selection, as is the case here. The effect size of mutations underlying traits under stabilizing selection within populations is expected to be smaller than the effect sizes of mutations in the early stages of an adaptive walk (Orr 1998; Remington 2015).

In addition to identifying regulatory variation at specific genes associated with body mass, we also used a systems biology approach to identify co-expression patterns associated with body size variation in wild mice. Gene co-expression networks capture biologically relevant relationships between genes that can be useful for understanding gene functions and interactions. Here, we have used this information to characterize co-expression modules that were associated with body size and identified regulatory variation within these co-expressed gene sets that may play a role in body size variation.

The tendency for body size to increase with latitude (i.e., Bergmann's rule) has been documented in many species, including humans (Ashton et al. 2000; Ruff 2002; Foster and Collard 2013), and reflects an evolved response to differences in temperature (Bergmann 1847). In humans, many candidate genes for metabolic disorders, such as obesity, also show evidence of climatic adaptation (Hancock et al. 2008). In humans, both ADAM17 and BCAT2 have been implicated in metabolic disease (Arribas and Esselens 2009; Newgard et al. 2009; Junyent et al. 2010; Menghini et al. 2013), and variation at ADAM17 has been identified in genome scans for selection (Pickrell et al. 2009; Parnell et al. 2010; Fumagalli et al. 2011) in addition to its association with body weight and obesity risk (Junyent et al. 2010).

Finally, this study provides evidence for the role of cis-regulatory variation in environmental adaptation in natural populations. While cis-regulatory variation has long been hypothesized to play a major role in adaptive phenotypic evolution, connecting regulatory variation with adaptive organismal phenotypes remains tricky. Combining eQTL mapping with genomic scans, as was done here, may be a fruitful approach for identifying adaptive regulatory variation in other natural systems.

Methods

Sampling

Mice used in this study were collected from five sampling locations (Supplemental Table S1; Supplemental File S1) along a latitudinal gradient in the eastern United States. Mice were sacrificed in the field and measurements (body weight, total body length, tail length) were taken at time of collection. Body mass index (BMI) was calculated as body weight/length2 (g/mm2). Liver tissue was collected in RNAlater and stored at 4°C overnight and then frozen to −80°C until RNA extraction with the Qiagen's RNeasy Mini kit.

mRNA-sequencing and mapping

For each sample, 100-base pair paired-end reads were sequenced on the Illumina HiSeq 4000 platform. RNA-seq reads were mapped with TopHat2 (Kim et al. 2013) to personal reference genomes, created by inserting variants into the mouse reference (GRCm38) and masking indels (see Supplemental Methods). We removed genes with fewer than 500 reads across samples (i.e., an average of 10 reads per sample). Gene expression was then quantile normalized and corrected for hidden factors and known covariates (individual sex and the first six principle components from genotype data to account for population structure) using a Bayesian approach (Supplemental Figs. S19, S20; Stegle et al. 2010, 2012).

Exome capture sequencing and identification of clinal outliers

The exome-sequence data was used to identify clinal outliers (Phifer-Rixey et al. 2018; see also Supplemental Methods). Libraries were enriched for exonic target regions and subsequently 100-bp paired-end reads were sequenced on the Illumina HiSeq 2000 platform, resulting in 2 GB of raw sequence data per individual. Forty-one of the 50 individuals for which there is exome-sequence data have matched RNA-seq libraries (see Supplemental Table S1). Reads were mapped with Bowtie 2 (Langmead and Salzberg 2012), and allele frequencies were estimated with ANGSD (Korneliussen et al. 2014). LFMM (Frichot et al. 2013) was used to identify covariance between environmental and genetic variation (see Supplemental Methods).

Cis-eQTL discovery

We performed cis-eQTL mapping using variant calls from RNA-seq and exome data (see Supplemental Methods). One limitation of this method is that the genotyping data set is limited to sites represented by these data (i.e., variant calls are largely limited to exomic regions of the genome). Consequently, many causal sites may not be typed and variants associated with expression may be tagging causal sites in LD. For the exome data set, depth per site of the targeted exome was ∼15×. For genes represented in the analysis, on average per individual we had sufficient coverage for ∼32% of bases within gene boundaries and ∼15% of bases in the 200-kb boundary used as the cut-off for cis-eQTL mapping (Supplemental Table S8).

To identify cis-acting eQTLs, we used a linear mixed model applied in the program GEMMA (Zhou and Stephens 2012) on expression residuals to associate expression with sequence variants (see Supplemental Methods). A relatedness matrix was computed and included as a covariate. We retained the variant with the lowest p-value for each gene and then performed a Bonferroni's correction. Variants with Bonferroni-corrected P-values of <0.05 were considered significant.

Weighted gene co-expression network analysis

We carried out a weighted gene co-expression network analysis on expression residuals following WGCNA protocols (Langfelder and Horvath 2008) to create expression modules. Each module is summarized by a representative eigengene, the first principle component of a given module. Each gene's expression was correlated with the module eigengene as a measure of the gene's centrality to the module, called module membership.

Data access

Illumina sequencing data from this study have been submitted to the NCBI BioProject (https://www.ncbi.nlm.nih.gov/bioproject) under accession number PRJNA407812. Museum accession numbers for samples used in this study are available in Supplemental File S1.