Introduction

Understanding the interactions between organisms and their environments is becoming increasingly important as environmental change and human activity change the original distribution of many species (Visser 2008; Wingfield et al. 2011; McCairns et al. 2016; Oostra et al. 2018). While some organisms can easily cope with environmental change, others run the risk of local or global extinction. An organismal response to environmental change partially depends on its ability of phenotypic change (Pfennig et al. 2010; Scoville & Pfrender 2010; Murren et al. 2015; Fox et al. 2019). For organisms that have recently colonized a new environment, phenotypic change may involve two stages. At the early stage, species may rapidly change their phenotype through plasticity without involving genetic changes. Under a persisting strong selective pressure, selection may bring genetic change to canalize phenotypic changes. Consequently, phenotypic plasticity has been considered to play an important stepping stone for genetically based evolutionary change at a late stage, so called “plasticity first” hypothesis (i.e., West-Eberhard 2003; Schwander & Leimar 2011; Levis & Pfennig 2016; Corl et al. 2018).

Under this hypothesis, plasticity may serve adaptive evolution in two alternative ways (Ho & Zhang 2018; Ho et al. 2020). On the one hand, if plasticity changes the phenotype in the same direction as adaptive evolution does, genetic variation can strengthen the change towards an optimal phenotypic value in the new environment (reinforcement plasticity). On the other hand, if the plastic change works in the opposite direction as that driven by adaptive evolution, the subsequent genetic change would have to revert the initial plastic response (reversion/reduction plasticity). As previous studies on morphological and physiological traits have revealed mixed evidence for the reinforcement and the reversion of plasticity (Pigliucci et al. 2006; Lande 2009; Pfennig et al. 2010; Moczek et al. 2011; Levis & Pfennig 2016; Corl et al. 2018), it remains unclear if and how genetic change acts to reinforce or reverse plasticity in the wild species (i.e., Campbell-Staton et al. 2021). The potential challenges are how to quantify phenotypic change relative to plasticity and adaptive evolution, and how to define what genes are involved in the phenotypic change, especially as many phenotypes are likely polygenically mediated (Novembre & Barton 2018; Pritchard et al. 2010).

Comparison of gene expression provides a potential approach to quantify the relative contributions of phenotypic plasticity and genetic variation since gene expression can bridge an organisms’ genotype to its cellular biology and, by extension, higher order physiological processes (Wray et al. 2003; Carrol 2008; Morris et al. 2014). While entire transcriptional program is orchestrated by the gene regulatory network, it is possible to trace the genes that are co-expressed or underpin the phenotypic change (Fukao et al. 2011). In addition, gene expression evolves in a stabilizing manner, where regulatory elements accumulate mutations that keep gene expression at an optimal level for physiological functions (Coolon et al. 2014; Gilad et al. 2006; Hodgins-Davis et al. 2015). Herein we used gene expression data to quantify regulatory plasticity and adaptation, as well as the relevant genetic changes in a group of wild birds that have recently colonized a highland environment, the Eurasian Tree Sparrows, Passer montanus (hereafter referred as tree sparrows).

The tree sparrow is a human commensal that has spread over a wide variety of habitats in the Eurasian continent, and that has been successfully introduced to Australia and North America (Widmann 1889; Graham et al. 2011). This species colonized the Qinghai-Tibet plateau approximately 2600 years ago (Qu et al. 2020), possibly concurrently with the introduction of barley agriculture in the highland. The Qinghai-Tibet plateau is a harsh environment with an average elevation of 4500 m above sea level (m a.s.l.) where only a few endemic animals can survive. For mammals and birds living at high elevations, hypoxia is one of the strongest selective pressures that drive physiological changes (Barve et al. 2016; McClelland & Scott, 2019). Compared to the lowland ancestral population, highland tree sparrows have evolved physiological and gene expression change for hypoxia tolerance (Sun et al. 2016; Qu et al. 2020), but how these changes relate to plasticity and adaptive evolution are largely unknown.

In order to explore this, we acclimated a group of lowland birds to an experimentally implemented hypoxic condition similar to that where we collected highland birds (i.e., 3200 m a.s.l). We could assess gene expression plasticity and regulatory adaptation in hypoxia tolerance by studying tree sparrows collected a) in their original lowland environment (representing the physiologically ancestral stage), b) shortly after their exposure to the new environmental stimulus (representing the ancestral plasticity at the initially plastic stage), and c) after having adapted to the highland environment (the colonized stage). Using the framework described in Figure 1, we quantified gene expression change stemming from hypoxia acclimation and highland colonization and defined genes with the expression plasticity being reinforced or reversed at the colonization stage. We test whether these genes increase genetic variation between ancestral and colonized populations using a permutation test. We use this framework in four datasets including groups of co-expressed genes associated to hypoxia tolerance and those correlated to muscle phenotypes for the cardiac and flight muscle, respectively. These independent comparisons congruently show that gene expression at the plastic stage is more often reversed than reinforced at the colonized stage, and the selection-driven genetic change depends on the magnitude of reinforcement and/or reversion plasticity.

Figure 1.

Results

Co-expressed genes involved in hypoxia tolerance generally reverse gene expression plasticity

We sequenced transcriptomic data of the cardiac and flight muscles from the lowland tree sparrows that had been kept under the hypoxia acclimation in this study (n=7), and re-analyzed transcriptomic data in the lowland and highland tree sparrows collected by Qu et al. (2020)(n=14, Supplementary Table 1). We focused on the flight and cardiac muscles because they affect thermogenic capacity and oxygen transport and delivery and are thus critical for the highland adaptation (Scott et al. 2015; Storz et al. 2010). Considering that response to environmental change often involves complex coordinated biological response, such as the co-regulation of genes underpinning physiological functions, we analyzed the resultant expression data for all individuals using weighted gene correlation network analysis (WGCNA, Langfelder & Horvath 2008). Weighted correlation network is a system biology method for identifying groups of highly co-expressed genes (modules), summarizing module-level expression and identifying important genes within modules (hub genes). Consequently, these co-expressed genes can capture expression-stage correlations due to the hypoxia-dependent expression. We identified six and nine modules of co-expressed genes with their expression changes that are significantly associated with the three stages (i.e., ancestral, plastic and colonized stages) for the flight and cardiac muscle, respectively (Figure 2a and Supplementary Figure S1). Across these modules, we identified 2413 and 508 hub genes based on the gene significance function (GS>1^st quartile and P<0.05) for the flight and cardiac muscles, respectively (Figure 2b-c).

Figure 2.

Using these candidate genes we examined how gene expression plasticity changed across the three stages. We compared the gene expression levels between the birds at the ancestral stage and those at the plastic stage, as well as between the birds at the ancestral stage and those at the colonized stage. We attributed genes showing the same direction of change to the reinforcement group and those showing opposite direction of change to the reversion group as described in Ho and Zhang (2018). Using a threshold of >50% increase/reduction of gene expression between two comparisons, we attributed 36 and 15 genes to the reinforcement group and 209 and 60 genes to the reversion group for flight and cardiac muscle, respectively. Interestingly, we found fractions of genes with reversion plasticity are larger than those of genes with reinforcement plasticity (Figure 2d). To confirm the robustness of this observation with respect to random sampling errors, we carried out a parametric bootstrap procedure as described in Ho and Zhang (2019), which aimed to identify genes resulting from genuine differences rather than random sampling errors. Bootstrap results also confirmed that genes exhibiting reversing plasticity significantly outnumber those exhibiting reinforcing plasticity (Figure 2d).

We then explore whether this result is robust with different magnitudes of plasticity by using a range of thresholds of <150% and >100%, <200% and >150%, and >200% to attribute gene to reinforcement and reversion groups. These different categories consistently showed an excess of genes with reversion plasticity (six out of eight comparisons, two-tailed binomial test, P<0.05, Fig. 2e-f). To further confirm the robustness of the discovered pattern against the arbitrariness of threshold selection, we also adopted different categorizations according the magnitude of plasticity (i.e., 20%, 40% and 60% bin settings along the spectrum of the reinforcement/reversion plasticity), and these analyses showed similar results (Supplementary Figure S2). Altogether, these results suggest that high-elevation colonization of tree sparrows generally reverses gene expression plasticity.

If the reinforcement and/or reversion gene expression plasticity in hypoxic tolerance were the target of selection in the high-elevation environment, we would expect to see an increase in the genetic divergence between the ancestral and colonized populations. To evaluate this we searched for genomic signatures of selection using resequencing data from 12 lowland and 11 highland tree sparrows generated from a previous study (Qu et al. 2020). Significantly elevated genetic divergence (i.e., F_ST) within the candidate genes beyond the genic background would provide evidence for selection acting on this group of co-expressed genes (i.e., polygenic adaptation). We thus calculated average F_ST for the genes with reinforcement or reversion plasticity using all SNPs included in the genic region, 2kb up-stream and down-stream regions, respectively. We considered the up-stream and down-stream regions because gene regulation may refer to genetic variation in nearby non-coding regions. By comparing to the F_ST distributions generated from 100 random samplings of the same number of SNPs from the genic background (see Methods), we found that the empirical F_ST values significantly increased in the genes with the levels of reinforcement/reversion plasticity reaching above 200% (P < 0.05, Figure 2e-f). In all cases, the empirical F_ST values increased significantly in the 2kb up-stream and/or down-stream regions of the candidate genes, but not in the genic regions. This observation thus suggests that selection-driven genetic change is more targeted on the non-coding regions and depends on the magnitude of reinforcement/reversion plasticity.

Gene expression and muscle phenotype analyses show adaptive plasticity and maladaptive plasticity

The evolutionary outcome of the phenotypic plasticity for a population colonizing a novel environment depends on the direction of plastic change with respect to the local optimum of novel environment. If the plastic change is close to the local optimum, i.e., adaptive plasticity, natural selection likely reinforces adaptive plasticity. Conversely, if the plasticity moves the phenotypic change away from the local optimum, i.e., maladaptive plasticity, natural selection should reduce or reverse the reaction norm and restore the phenotypic change back to the original ancestral values (Campbell-Staton et al. 2021). In order to investigate phenotypic plasticity and its evolutionary consequence, we analyzed two datasets that included gene expression and muscle phenotypes collected from the lowland ancestral and highland colonized tree sparrows (Qu et al. 2020).

Based on the correlation analysis between gene expression and muscle phenotypic values of the ancestral and colonized tree sparrows, we categorized the direction of expression-phenotype correlation as positive or negative regulators for the muscle phenotype-associated genes (P < 0.05, Figure 3a). We identified 2037 and 1866 muscle phenotype-associated genes in the flight and cardiac muscle, respectively (Figure 3b-c). Out of these genes we found 328 positive regulators and 1709 negative regulators for the flight muscle, and 843 positive regulators and 1023 negative regulators for cardiac muscle (Supplementary Figure S3). We then characterized the differences in gene expression observed for the ancestral-plastic groups and the ancestral-colonized groups as congruent or incongruent with increased hypoxia tolerance (see Methods). Because the muscle phenotypic changes observed in the colonized tree sparrows have been shown to enhance oxygen delivery and improve metabolic capacity in the highland animals (Giordano 2005; Storz et al. 2010; Scott 2011; Scott et al. 2015), we deemed that the gene expressions in which the direction of plastic change matched the expectation for increase hypoxia tolerance as putatively adaptive, while those that opposed this expectation were considered to be putatively maladaptive (Figure 3a). Using a threshold of >50% of plasticity change, we found that the tree sparrows displayed more genes with maladaptive plasticity than genes with adaptive plasticity (two-tailed binomial test, flight muscle, 198 vs. 48, P<2.2e^-16, cardiac muscle 186 vs. 166, P<0.05). These results are robust with random sampling errors as our parametric bootstrap analyses also revealed similar results (Figure 3d).

Figure 3.

We subsequently examined expression-stage relationship using the thresholds (<150% and >100%, <200% and >150%, and >200%) and bin settings (20%, 40% and 60% bin settings) as above-mentioned to explore the magnitudes of adaptive and maladaptive plasticity. We found that the fractions of genes with maladaptive plasticity were larger than those of genes with adaptive plasticity (two tailed binomial test, P<0.01, Figure 3d). This observation is robust with different thresholds and bin settings for the flight muscle, but slightly sensitive to bin setting for the cardiac muscle, indicating by more than 50% of comparisons supporting an excess of genes with maladaptive plasticity (Supplementary Figure S4). These results suggest that gene expression plasticity induced by hypoxia exposure at the early stage may be mainly maladaptive for muscle physiological performance in highland environment.

The observed adaptive and maladaptive gene expressions in the tree sparrows suggest that selection may have acted to reinforce and/or reverse gene expression plasticity. We then quantified the genetic divergence of these genes between ancestral and colonized tree sparrows by testing for statistically significant increase of F_ST as compared to the F_ST distributions generated from 100 random samplings (see Methods). We found that the empirical F_ST values were significantly larger in the genes within the <200% and >150% category (cardiac muscle, P<0.05) and >200% category (flight and cardiac muscle, P<0.05), but not in those within other categories with low levels of reinforcement and reversion plasticity. Additionally, we found that genetic divergence between the ancestral and colonized tree sparrows tends to increase in the 2kb up-stream and down-stream regions of these genes (P < 0.05, Figure 3e-f). These results suggest that selection on regulatory adaptation (i.e., genetic divergence) may only be expected if large strength of selection is needed to bring initial plasticity toward adaptive optimum of the novel environment.

Discussion

An understanding of phenotypic plasticity and its consequence for adaptive evolution is important in evolutionary biology but remains challenging (Ghalambor et al. 2007; 2015; Campbell-Staton et al. 2021; Kenkel & Matz 2017; Rivera et al. 2021). Using an integrative approach of field-based studies, experimentally implemented acclimation experiment, multiomic and muscle phenotypic data, our study demonstrates the importance of gene expression plasticity in facilitating initial survival and population persistence at early stage of colonization. In particular, our results provide novel insight into an issue of active debate within evolutionary biology (i.e., reinforcement vs. reversion plasticity) by showing that the evolutionary consequence of plasticity depends on the adaptive plasticity (reinforcement) or maladaptive plasticity (reversion). We also demonstrate that natural selection drives genetic divergence when the magnitude of reinforcement and/or reversion plasticity becomes intense, providing novel insights into the mechanism of plasticity first and gene as followers (Schwander & Leimar 2011).

Our gene expression and sequence analyses provide several lines of evidence in support for gene expression plasticity contribution to the high-elevation colonization of tree sparrows. First, we show that co-expressed genes and muscle phenotype associated genes change either the intensity or the direction of their expression profiles, and that gene expressions at the plastic stage are more often reversed than reinforced at the colonization stage. Second, by combining gene expression and muscle phenotype, we show that selection operates intensely on maladaptive plasticity during the initial stage of responding novel environment (Ghalambor et al. 2007; 2015; Ho and Zhang 2018; Campbell-Staton et al. 2021). Third, the genes with intensified levels of reinforcement and reversion expression plasticity show large genetic divergence (F_ST) in their up-stream and down-stream regions between the ancestral and colonized populations, suggesting that selection has driven genetic changes in the noncoding regions after the tree sparrows had colonized the novel environment.

These observations thus support that plasticity serves as a stepping stone in a successful colonization of a new environment and that genetic changes that evolve at a late stage will modify the plastic change to an optimal phenotype in the new environment (Corl et al. 2018; Schwander & Leimar, 2011). It is interesting to think about why phenotypic change at the plastic stage sometimes differs from those at the colonized stage. A likely reason is that initially plastic change allows organisms to survive upon a sudden environmental shift but the fitness is much reduced compared to that after a long-term colonization to the novel environment (Fischer et al. 2016; Huang & Agrawal 2016; Leonard & Lancaster 2020; Kuo et al. 2023). Thus, the overall physiological states of the organisms right after an environmental shift may result in low fitness. To coordinate the phenotype values to close the organismal fitness, a genetic change may be required modifying plastic change of phenotype to approach the optimum to the novel environment (Storz 2021). Thus, the consequence of genetic change on plasticity is dependent on the cost of plastic change (Ho and Zhang 2018; Corl et al. 2018).

A limitation of our study is that we could only consider hypoxic condition in the acclimation experiment, while also other factors associated with high-elevation environments, e.g., low temperature and UV radiation, may be of importance. Although we assume that hypoxia is the dominant selective force in the Qinghai-Tibet Plateau (see also Sun et al. 2016; Qu et al. 2020), the use of a single environmental stimulus makes our conclusions more conservative. Also, it is reasonable to assume that plastic change in organisms is time-dependent and that our results cannot catch gene expression occurring at later stage of plastic response (i.e., after one month). Since we could only run the acclimation experiment for one month it is likely we have not identified all genes related to hypoxia tolerance. Despite this, we could observe a clear pattern of reinforcement and reversion plasticity and the associated genetic changes. Our sample sizes are rather small (n=12 for the flight muscle and n=9 for the cardiac muscle), mostly because of the logistical challenge of keeping tree sparrows in the necessary common garden experiments. In order to compensate for this, we utilized comparative transcriptomics to generate the four independent datasets (a group of co-expressed genes and a group of phenotype-associated genes for the flight and cardiac muscles, respectively). As the analyses of these four independent datasets show similar results, we believe the conclusions drawn from our study are robust.

Despite some limitations, our study demonstrates an easy implemented framework for quantifying phenotype plasticity and how to relate this to genetic change. Short-term acclimation response of lowland populations has been found to be different from the populations adapted to the high-elevation environments. For example, people newly exposed to high-elevation environment increase their red blood cell as compared to people at sea level, which is likely mediated by HIF and EPO genes in hypoxia response cascade (Erzurum et al. 2007; Peng et al. 2017). This physiological change, however, is less activated in native people at high-elevation when compared to people at sea level, demonstrating an example of reversion of phenotypic plasticity (Storz 2021). Together with previous studies, our work shows that short-term acclimation response and long-term adaptive evolution in hypoxia response may affect different evolutionary pathways. Our work provides understanding how plasticity facilitates species to invade new habitats and survive environmental change. Such an understanding is important, in particular when numerous species are being introduced to new regions of the globe through human intervention and spread as invasive species (Visser 2008).

Materials and Methods

Sampling

We generated transcriptomic data from hypoxia-acclimated lowland tree sparrows by exposing five lowland tree sparrows to hypoxic condition for 30 days using a hypoxic chamber (Qu et al. 2020). The oxygen content in the chamber was set to 14% of the oxygen content (simulating the oxygen concentration at 3200 m a.s.l., 70% of ∼20%, latter of which is the content of oxygen content at sea level), which is similar to condition where we collected the highland tree sparrows (Qinghai Lake, 3200 m a.s.l.). We collected flight muscle from the four individuals and cardiac muscle from three individuals. In addition, we extended the analyses of the transcriptomic and muscle histological data from tree sparrows that were collected from lowland (Beijing and Hebei, 100 m a.s.l.) and highland (Qinghai Lake, 3200 m a.s.l., Figure 1) collected in previous study (Qu et al. 2020). In total, we used twelve samples for flight muscle and nine samples for cardiac muscle (Supplementary Table 1). The ancestral lowland and hypoxia-exposed lowland tree sparrows were collected from the same locality (Beijing) and the same season (pre-breeding). Only adult birds with similar body weight (approximately 18g) were used.

Transcriptome sequencing

RNA libraries were constructed for the flight and cardiac muscles and sequenced on an Illumina HiSeq4000 platform. After filtering low-quality, adapter-contaminated, and N-rich reads (>10%), a total of 134, 128 and 120 million reads of transcriptional data were generated for the lowland, hypoxia-exposed lowland tree sparrows and highland tree sparrows, respectively (Supplementary Table S1). We mapped cleaned reads against the tree sparrow genome using STAR (Dobin et al. 2013). We compared log-transformed transcript per million (TPM) across all genes and determined the most conserved genes (i.e., coefficient of variance ≤ 0.3 and average TPM ≥ 1 for each sample) for the flight and cardiac muscles, respectively (Hao et al. 2023). We then compared the median expression levels of these conserved genes and found no difference among the lowland, hypoxia-exposed lowland, and highland tree sparrows (Wilcoxon signed-rank test, P<0.05, Supplementary Figure S5), suggesting that batch effect had little influence on the transcriptomic data. We then used TPM values to calculate gene expression level and intensity using RSEM (Li & Dewey 2011).

WGCNA analyses

We used WGCNA v. 1.41-1 (Langfelder and Horvath, 2008) to identify regulatory architecture for the flight muscle and cardiac muscle transcriptomes. Specifically, we used a principal component analysis (PCA) to summarize modules of genes expression with blockwiseModules function, and then used module eigengene values of the first principal component (PC) to test correlation between module expression and stages. We used Pearson linear correlation (cor function) to identify genes with expression reinforcement and nonlinear regression (polynomial regression model: Model = lm(y ∼ poly(x, 5, raw = TRUE), data=data) to identify genes with expression reversion. Regulatory modules were identified as branches of the resulting cluster tree via the dynamic tree-cutting method and highly correlated modules (P<0.1) were merged. We used “GS > 1^st quartile of GS” and “P<0.05” to identify important genes within modules (hub genes).

Identify gene regulation associated with muscle phenotype

To identify the genes with expression levels that correlated to evolved differences in the muscle phenotypes in the tree sparrows, we used the lowland and highland individuals from Qu et al. (2020) for which RNAseq and muscle phenotypic data were available. We correlated the levels of gene expression with muscle phenotypes in the lowland and highland individuals. We used three phenotypes for the flight muscle (capillary per fiber, fiber area and perimeter) and two for the cardiac muscle (capillary density and fiber diameter). We used Pearson correlation to test the association between each of the phenotypic traits and expression profile of each gene (TPM). A threshold of P<0.05 was used to detect muscle-associated genes. We then followed Campbell-Staton et al. (2021) to categorize genes as positive (positively correlated with muscle phenotypes) or negative regulators (negatively correlated with muscle phenotypes).

Subsequently, we compared the direction of gene expression plasticity at the plastic stage with those of the positive and negative regulators. Briefly, the direction of gene expression plasticity at the plastic stage matched the expectation for positive regulator (increased expression) or negative regulator (decreased expression) in the colonized tree sparrows was considered to indicate adaptive plasticity. Conversely, those that opposed this expectation were considered to be maladaptive, i.e., the direction of gene expression plasticity at the plastic stage showed decreased expression for the positive regulators and increased expression for negative regulators (Figure 3a).

Reinforcement and reversion plasticity analyses

To investigate the magnitude of reinforcement and reversion plasticity, we followed the method described in Ho and Zhang (2018). Specifically, expression levels of each gene in the three stages were treated as E_lowland (expression level at ancestral stage), E_hypoxia (expression level at plastic stage) and E_highland (expression level at colonized stage). We used a threshold of 50% of the ancestral gene expression level, i.e., 50% * E_lowland, to detect an excess of gene expression change. We identified genes with an excess of plastic change if it satisfied the condition of |(E_hypoxia-E_lowland)|> 50% * E_lowland. Likewise, we identified genes with an excess of evolutionary change if it met the condition of |(E_highland-E_hypoxia)|> 50%* E_lowland. If the genes showed same direction in the expression change at the plastic and colonized stages, we regarded these genes as reinforcement expression. If the genes showed opposite direction in the expression change at the plastic and colonized stages, we regarded these genes as reversion expression. We categorized these genes to reinforcement and reversion groups and then tested if the proportions of the genes in the two groups were significantly different from the expected proportions using two-tailed binomial tests.

To explore whether an excess of genes with reversion plasticity was subject to random sampling errors (Mallard et al. 2018), we used a parametric bootstrap method as described in Ho and Zhang (2019). Specifically, we simulated the mean expression level of a gene at each stage following a Gaussian distribution within the mean equal to the observed mean expression of this gene at this stage and the standard deviation equal to the estimated standard error. We then draw a random variable from the above Gaussian distribution to represent an observation of the mean expression level of this gene at this stage. We drew random variables representing the mean expression level of gene at each of the three stages (i.e., ancestral, plastic and colonized stages), and then computed gene expression changes at plastic stage and colonized stage and determined the reinforcement or reversion plasticity. This process is repeated 1000 times for each gene. If this gene exhibited reinforcement or reversion plasticity at least 950 repeats (P<0.05), we considered this gene showing reinforcement or reversion plasticity. The Gaussian distribution and random sampling were generated using rnorm function of R.

To test the impact of the choice of thresholds on the results, we analyzed the data using a range of thresholds, i.e., <100% and >50%, <150% and >100%, <200% and >150%, and >200%. This setting aimed to obtain independent genesets for subsequently genetic divergence analyses, because a geneset selected by relax threshold (i.e., 50%) also included genesets obtained by the stringent threshold (i.e., 100%, 150% and 200%). In addition, we also implemented other categorization scheme to check the robustness of the results.

Specifically, we used three different bin settings (i.e., 20%, 40% and 60% bin settings along the spectrum of the reinforcement/reversion plasticity), to group genes according their magnitude of plasticity. For each category in the three bin settings, we compared whether the proportion of the genes with reversion plasticity differed from that with reinforcement plasticity using two-tailed binomial tests.

Comparing genetic divergence for genes with reinforcement and reversion expression plasticity

To calculate the genetic divergence between the lowland ancestral and highland colonized tree sparrows, we used 12 lowland and 11 highland individuals for which re-sequencing data were available (Qu et al. 2020). After mapping the raw reads to the tree sparrow genome using BWA v0.7.17 (Li & Durbin 2009), we obtained mean sequence coverage of 17× for each individual (Table S2). We called variants with GATK v 3.7 (McKenna et al. 2010) and Samtools v1.2 (http://www.htskib.org/) and filtered single nucleotide polymorphism (SNP) using VCF tools and GATK with minimum coverage = 138, root mean square mapping quality >= 20, distance of adjacent SNPs >=5, distance to a gap >=5bp and read quality value >=30.

We calculated average F_ST using Vcftools (Danecek et al. 2011) with SNPs included in the genic region, 2kb up-stream and down-stream regions of the genes within each of the four categories, i.e., <100% and >50%, <150% and >100%, <200% and >150%, and >200%. The up-stream and down-stream regions were considered because gene regulation may refer to genetic variation in the nearby genic regions. We permuted F_ST distributions for the genes within each category by random sampling the same number of SNPs, allowing for a fluctuation of 5% of total SNPs. We compared the empirical F_ST values to permutated F_ST distributions generated from 100 samplings and considered a threshold of P<0.05 to be statistically significant (empirical F_ST>95% percentile of permutated F_ST distribution). We calculated empirical F_ST values and permutated F_ST distributions for the genic region, 2kb up-stream and down-stream regions separately. To avoid the sampling bias (i.e., sampling SNPs from neutral regions), our SNP sampling was constrained to genic region, 2k up-stream and down-stream regions of all genes (i.e., 16925 genes). All analyses were conducted with the R statistical software package (R Foundation for Statistical Computing 2017).

Data accessibility

DNA sequencing reads used in this study have been deposited in Short Read Archive under the project number PRJNA417520.

Authors’ contributions

Q.Y. designed research; S.H., Y.H., X.L., S.G. and Q.Y. performed research; F.L. provided highland samples; Q.Y. and S.H. wrote the paper.

Acknowledgements

We acknowledge Ying Xiong for logistic work in the hypoxia-exposed experiment. This research was funded by the Third Xinjiang Scientific Expedition and Research (XIKK) (2022xjkk0205), Second Tibetan Plateau Scientific Expedition and Research (2019QZKK0501) and the National Natural Science Foundation of China (NSFC32020103005).