Phenotypic plasticity co-varies with elevational range in two avian species of elevational migrants in the Himalayas

Xue, Boning; She, Huishang; Zuo, Lirong; Bates, John M.; Cheviron, Zachary A.; Wu, Yongjie; Li, Dongming; DuBay, Shane G.; Qu, Yanhua

doi:10.1038/s41467-025-60770-w

Download PDF

Article
Open access
Published: 18 June 2025

Phenotypic plasticity co-varies with elevational range in two avian species of elevational migrants in the Himalayas

Nature Communications volume 16, Article number: 5316 (2025) Cite this article

2422 Accesses
2 Altmetric
Metrics details

Subjects

Abstract

Montane birds represent a large proportion of global avian biodiversity and are highly sensitive to seasonal climatic shifts. Many montane species undertake elevational migration in response to seasonal climate change, but how they cope with disparate environmental pressures along the elevational gradient is poorly understood. Here, we investigate phenotypic plasticity in two closely related elevational migrants that differ in their seasonal movements (small-scale vs. large-scale elevational migrants). Through common-garden experiments with adult male birds of these species, we compared gene and protein expression plasticity in response to changes in temperature and partial pressure of oxygen (PO₂). The small-scale elevational migrant exhibited greater plasticity in response to changes in temperature, and the large-scale elevational migrant exhibited greater plasticity in response to changes in PO₂. These findings correspond to annual variation in environmental conditions experienced by each species, suggesting that phenotypic plasticity may co-vary with the seasonal elevational distributions of these elevational migrants.

Adaptation to climate change through dispersal and inherited timing in an avian migrant

Article 14 September 2023

Pursuit and escape drive fine-scale movement variation during migration in a temperate alpine ungulate

Article Open access 02 July 2024

The influence of climate variability on demographic rates of avian Afro-palearctic migrants

Article Open access 16 October 2020

Introduction

Twice a year, millions of birds migrate across varying distances, moving between breeding and nonbreeding grounds; avian migration is an extraordinary behavior that enables individuals to seasonally shift environments^1,2. During these largely latitudinal movements, migrants take days and even months to adjust to changes in day length, temperature, and other environmental factors that they face on flyways^3,4,5. Unlike latitudinal migrants, montane birds often migrate seasonally across elevation, moving up and down elevational gradients, a process known as elevational (or altitudinal) migration^6,7. Despite sometimes only travelling short geographic distances, elevational migrants can face dramatic changes in temperature and partial pressure of oxygen (PO₂) in a span of hours to a few days. For example, montane birds that undergo elevational shifts of ≥2000 m between seasons have to cope with changes in PO₂ by ≥20%⁸. Given the physiological challenges associated with shifts in elevation, elevational migrants exhibit phenotypic plasticity to facilitate acclimation to PO₂ and temperature variation experienced during the annual cycle⁹. For example, the Southern Giant Hummingbird (Patagona gigas) exhibits striking plasticity in blood traits related to seasonal elevational shifts, such as in hemoglobin concentration and hematocrit¹⁰, and the Bar-headed Goose (Anser indicus) shows plastic increases in O₂ binding capacity that facilitates its annual movements between high and low elevations¹¹. However, despite these case studies, how phenotypic plasticity co-varies with seasonal elevational distributions of elevational migrants remains largely unknown.

Theory predicts that phenotypic plasticity should correlate with temporal variation in environmental selection pressures (e.g., temperature and PO₂ for montane species)^12,13. Plasticity enables an organism to optimize fitness across different environmental conditions^14,15,16. For example, montane species with narrow elevational ranges may require greater phenotypic plasticity to accommodate seasonal fluctuations in temperature, but less plasticity in response to PO₂ changes since they experience little variation in PO₂ throughout the year. In contrast, species with large elevational ranges may require greater phenotypic plasticity to cope with change in PO₂ across the elevational gradient, but less plasticity in response to change in temperature, as species undergoing the largest elevational changes often experience the least temperature variability¹⁷. Additionally, the seasonal fluctuations in temperature increase for species with ranges at higher elevations¹⁸. While a strong theoretical framework exists for understanding phenotypic plasticity in relation to seasonal elevational movements, empirical tests remain scarce.

To understand how seasonal elevational distributions influence phenotypic plasticity we investigate how plasticity co-varies with variation in temperature and PO₂ in two closely related bird species from the Sino-Himalayan mountain region, a biodiversity hotspot in Asia¹⁹. Avian habitats in the Himalayas span an elevation gradient of up to 5000 m above sea level (a.s.l.), and nearly 65% of high-elevation breeding birds in this region undergo elevational migration^20,21. One of our focal species, the Golden Bush-robin (Tarsiger chrysaeus), is a large-scale elevational migrant that breeds at elevations between approximately 3000 and 4600 m a.s.l. and spends non-breeding periods at elevations as low as 500 m a.s.l.²². Its closest relative in this region, the White-browed Bush-robin (Tarsiger indicus), is a small-scale elevational migrant that breeds at high elevation in sympatry with T. chrysaeus but spends non-breeding periods at much higher elevations, above 2000 m a.s.l. (Figs. 1 and 2a, the degree of shifts for the two species are based on previous studies on occurrence data of elevational migrants, e.g., Ref. ^22,23). Across seasons, T. chrysaeus experiences dramatic changes in PO₂ between breeding and nonbreeding elevations but relatively subtle temperature variation compared to T. indicus²⁴. In contrast, T. indicus experiences a relatively subtle change in PO₂ between breeding and nonbreeding elevations but a large seasonal shift in temperature at upper elevations compared to T. chrysaeus²⁴. These species-specific differences in seasonal elevational distributions provide a natural framework to test the relative and combined roles of temperature and PO₂ in shaping phenotypic plasticity.

**Fig. 1: A schematic of hypothetical plasticity of the two elevational migrants in response to changes in temperature and PO₂.**

**Fig. 2: An overview of gene expression plasticity between treatments and between species.**

We conducted common-garden experiments in which we manipulated temperature and PO₂ for both small-scale and large-scale elevational migrants, quantifying gene and protein expression differences in functionally important tissues for thermoregulation, respiration, and metabolism (i.e., the heart, flight muscle, lung and liver; a schematic framework is shown in Fig. 1). Gene and protein expression differences are of particular interest because phenotypic plasticity stems from environmentally sensitive gene expression and regulation of gene products²⁵. Acclimation to environmental change often involves complex coordinated biological responses with the entire transcriptional program orchestrated by gene regulatory networks. It is thus possible to trace phenotypic changes that are co-expressed at optimal levels for physiological functions, gaining a comprehensive and quantitative understanding of plasticity^26,27.

We predict that the seasonal elevational ranges of the two species shape their capacity for phenotypic plasticity. Specifically, the small-scale elevational migrant experiences greater seasonal variation in temperature and less seasonal PO₂ variation compared to the large-scale elevational migrant. Thus, we predict that the small-scale elevational migrant (T. indicus) should exhibit greater phenotypic plasticity (i.e., gene or protein expression difference) in response to changes in temperature but less plasticity in response to changes in PO₂. In contrast, we predict that the large-scale elevational migrant (T. chysaeus) will show greater plasticity in response to changes in PO₂ but less plasticity in response to change in temperature. Alternatively, because the small-scale elevational migrant experiences a wider range of temperatures, it may be robust to greater temperature regimes and might therefore exhibit a blunted gene expression response to temperature; whereas the large-scale elevational migrant may exhibit diminished plasticity (relative to the small-scale migrant) when exposed to dramatically different PO₂ conditions. Our results demonstrate that the small-scale elevational migrant generally exhibits a greater plastic response to change in temperature but less response to change in PO₂ than the large-scale elevational migrants. We find support for our prediction that phenotypic plasticity co-varies with the scale of environmental variation experienced by the species across the year, providing empirical support for long-standing theoretical predictions in the field.

Results

An overview of gene expression plasticity between treatments and between species

We acclimated the small-scale elevational migrant (T. indicus, n = 21) and the large-scale elevational migrant (T. chrysaeus, n = 23) to one of four experimental treatments, manipulating temperature and PO₂ (see section “Methods”, Fig. 2b). The treatments included: (1) cold/hypoxia (8 ± 2 °C and 92 mmHg), (2) warm/hypoxia (27 ± 2°C and 92 mmHg), (3) cold/normoxia (8 ± 2°C and 151 mmHg), and (4) warm/normoxia (27 ± 2 °C and 151 mmHg). Warm treatments correspond to a general thermoneutral zone, where adult birds maintain body temperature with minimal energy expenditure, while cold treatments represent temperatures below this zone, requiring increased metabolic activity for thermoregulation²⁸. Normoxia and hypoxia treatments simulate PO₂ conditions at approximately 500 and 4650 m a.s.l., respectively, approximating the species’ lower and upper elevations limits. After a two-week acclimation period, we collected heart, flight muscle, liver, and lung tissues for RNA-sequencing analyses (Supplementary Tables S1–S2 and Supplementary Data 1).

In both species, we first analyzed transcriptomic data to test for genetic variation between individuals across treatments, as such variation could confound results by influencing gene expression differences. We found no significant differences in nucleotide diversity or genetic heterozygosity among treatments or sampling years (Mann–Whitney U test, all comparisons, P > 0.05, Supplementary Fig. S1). Given the absence of genetic differences between treatments, we then evaluated gene expression plasticity by identifying differentially expressed genes (DEGs) within each species. Specifically, we examined DEGs across three comparisons: (1) cold/hypoxia vs. warm/normoxia to assess the combined effects of temperature and PO₂ change, (2) cold/hypoxia vs. warm/hypoxia and cold/normoxia vs. warm/normoxia to assess the effects of temperature change alone, and (3) cold/hypoxia vs. cold/normoxia and warm/hypoxia vs. warm/normoxia to assess the effects of PO₂ change alone. In both species, temperature-PO₂ change and temperature change alone resulted in more DEGs than PO₂ change alone (Fig. 2c and Supplementary Table S3), suggesting that temperature change generated greater gene expression plasticity than PO₂ change.

We then compared gene expression plasticity between the two species. In response to combined change in temperature and PO₂, the small-scale elevational migrant exhibited a greater number of DEGs in the flight muscle and liver, whereas the large-scale elevational migrant showed more DEGs in the heart and lung tissues (two-tailed exact binomial test, P < 0.05, Fig. 2d). When assessing the effect of changes in temperature alone (cold/hypoxia vs. warm/hypoxia), the small-scale elevational migrant exhibited more DEGs than the large-scale elevational migrant in all tissues except the heart (P < 0.001, Fig. 2d). A similar pattern was observed in the cold/normoxia vs. warm/normoxia comparison, where the small-scale elevational migrant displayed greater number of DEGs in the heart, flight muscle and lung (P < 0.001, Fig. 2d). These findings indicate that the small-scale elevational migrant exhibits greater gene expression plasticity in response to temperature change, though this effect appears to be tissue-specific.

When analyzing the effect of PO₂ change alone (cold/hypoxia vs. cold/normoxia), we found that the large-scale elevational migrant exhibited more DEGs than the small-scale elevational migrant in the flight muscle and lung tissues (P < 0.001). Similarly, a parallel comparison between warm/hypoxia vs. warm/normoxia revealed the same pattern in the heart, liver, and lung tissues (P < 0.001, Fig. 2d). Taken together, these results suggest that the large-scale elevational migrant exhibits greater gene expression plasticity in response to PO₂ change than the small-scale elevational migrant.

Regardless of changes in temperature or PO₂, or their combined effects, the two species recruited different DEGs, as indicated by 0% to 7.7% shared DEGs between them (Supplementary Fig. S2). However, these genes were often enriched for similar biological functions. For example, in flight muscle and liver tissues, DEGs identified in response to temperature-PO₂ change and temperature change alone were enriched for gene ontology (GO) terms associated with metabolic activities, biological regulation and response to stimuli (Supplementary Data 2). Furthermore, DEGs identified in both species were frequently involved in the same pathways, such as glycerolipid metabolism, fatty acid metabolism, glycine, serine and threonine metabolism and N-Glycan biosynthesis (false discovery rate (FDR) adjusted P value < 0.1, Supplementary Fig. S3). These results suggest that, despite recruiting different genes, the two species may invoke similar functional responses within the same pathways in response to temperature and PO₂ changes.

Gene expression plasticity is greater in response to temperature change than PO₂ change

Gene expression plasticity in response to temperature and PO₂ changes may involve coordinated biological responses, such as the co-regulation of genes underpinning physiological functions. To explore this, we used weighted gene correlation network analysis (WGCNA)²⁹ to examine gene expression associated with the different treatments. To compare gene expression plasticity in response to temperature and PO₂ changes (combined and individually), we used all four treatments to identify groups of highly co-expressed genes (modules) and key genes within those modules (hub genes)³⁰. We identified ten co-expression modules in the heart, 16 in the flight muscle, 19 in the liver and 24 in the lung (Pearson correlation, P < 0.05, Supplementary Fig. S4). Within these modules, we detected 664 hub genes in the heart, 2568 in the flight muscle, 1638 in the liver and 676 in the lung, based on gene significance (Pearson correlation, P < 0.05, Supplementary Fig. S4).

Using the hub genes, we performed a vector analysis to summarize gene expression differences in a multivariate trait space, as described in refs. ^31,32. Gene expression divergence between treatments was quantified by vectors connecting the multivariate means (centroids) in pairwise comparisons (see Methods). Using the vectors from the cold/hypoxia treatment as the baseline (i.e., the environmental conditions corresponding to the breeding elevations of the two elevational migrants), we estimated vector lengths for warm/normoxia (combined changes in temperature and PO₂), warm/hypoxia (changes in temperature), and cold/normoxia (changes in PO₂). The estimated vector lengths were larger for warm/normoxia and warm/hypoxia compared to cold/normoxia (Fig. 3a and Supplementary Fig. S5). This pattern was consistent across all four tissues and both species (Fig. 3a and Supplementary Table S4).

**Fig. 3: Quantifying gene expression plasticity in the small-scale and large-scale elevational migrants using vector analysis.**

We then used the vector angle to estimate the direction of gene expression differences resulting from changes in temperature and PO₂ independently (Fig. 3b, see section “Methods”). We found that the angles between the cold/normoxia and warm/normoxia vectors deviated more from the baseline (i.e., cold/hypoxia) than those between the vectors of the warm/hypoxia and warm/normoxia (Fig. 3b and Supplementary Table S5). These results again suggest that gene expression plasticity in response to temperature change is greater than that produced from PO₂ change.

Species differ in gene expression plasticity in response to changes in temperature and PO₂

We next compared gene expression plasticity in response to changes in temperature and PO₂ individually between the two species. Using cold/hypoxia as the baseline, we first identified hub genes associated with warm/hypoxia to assess the effects of temperature change (Supplementary Fig. S6). This analysis yielded 156 hub genes in the heart, 1214 in the flight muscle, 983 in the liver, and 545 in the lungs. We summarized the expression changes of these hub genes using vector analysis and found that the vector lengths were greater for the flight muscle, liver and lung tissues in the small-scale elevational migrant than in the large-scale elevational migrant, but smaller for the heart (Mann–Whitney U test, P < 0.001, Fig. 3c and Supplementary Table S6).

We then compared gene expression plasticity in response to PO₂ change between species by identifying hub genes associated with cold/hypoxia and cold/normoxia (Supplementary Fig. S7). The hub genes identified for changes in PO₂ exhibit minimal overlap with those detected for changes in temperature (i.e., 3.2% to 15.7% hub genes were identified in both comparisons, Supplementary Fig. S8). Vector lengths between cold/hypoxia and cold/normoxia vectors were greater in the large-scale elevational migrant than in the small-scale elevational migrant in all organs except the liver (Mann–Whitney U test, P < 0.01, Fig. 2d and Supplementary Table S7). Of the four tissues, liver tissue had a smaller sample size than the minimum required for WGCNA analysis (n = 13 as compared to >15 required by WGCNA³³), so we used a jackknife resampling approach to test the robustness of liver analysis. We found that the jackknife resampling modules and full data modules showed considerable similarity (Dice similarity coefficient, 0.84 for cold/hypoxia vs. warm/hypoxia and 0.75 for cold/hypoxia vs. cold/normoxia; Jaccard index: 0.72 for cold/hypoxia vs. warm/hypoxia and 0.6 for cold/hypoxia vs. cold/normoxia, Supplementary Fig. S9). Additionally, vector analysis using jackknife resampling modules yielded similar results to those from the full modules (Supplementary Fig. S10), indicating that our findings are robust despite the modest sample size of the liver tissue. Collectively, these results suggest that the small-scale and large-scale elevational migrants exhibit differences in gene expression plasticity in response to changes in temperature and PO₂ independently.

Small-scale elevational migrant and large-scale elevational migrant exhibit plasticity differences in flight muscle glycolytic activity

To further examine plasticity between the two elevational migrants, we focused on genes related to metabolic activities, as most GO terms enriched for hub genes were related to these functions (average 58%, Supplementary Fig. S11 and Supplementary Data 3). Specifically, for each species we calculated the coefficient of variation (CV, an indicator of phenotypic plasticity)³⁴ of gene expression across the four acclimation treatments for hub genes that are functionally relevant to fatty acid oxidation (i.e., fatty acid metabolism pathway), aerobic metabolism (i.e., citrate cycle pathway), and glycolytic metabolism (i.e., glycolysis/gluconeogenesis pathway). In all four tissues, we found no species difference in CVs of hub genes associated with fatty acid oxidation or the citrate cycle. However, for glycolysis/gluconeogenesis, we observed that the small-scale elevational migrant exhibited a larger CV of gene expression in the flight muscle compared to the large-scale elevational migrant, but not in the other tissues (Fig. 4a and Supplementary Fig. S12). These results suggest that the small-scale elevational migrant has maintained a greater capacity for plasticity in flight muscle glycolysis than the large-scale elevational migrant.

**Fig. 4: Small-scale elevational migrant and large-scale elevational migrant show plasticity differences in flight muscle glycolytic activity.**

To further investigate changes in glycolytic enzyme activity in the flight muscle of the two species, we conducted a flight muscle proteomic analysis. We assessed protein expression plasticity across the four acclimation treatments by calculating CV for each protein in the glycolysis/gluconeogenesis pathway for the two species. Of the 14 proteins identified as hub genes in this pathway, we found that nine exhibited larger CVs in the small-scale migrant than in the large-scale migrant (Fig. 4b). Notably, two rate-limiting glycolytic enzymes, hexokinase (HK) and pyruvate kinase (PKM), displayed greater protein expression plasticity in the small-scale than in the large-scale migrant, indicating a higher capacity for modulating glycolytic activity. Together, these results suggest that the small-scale elevational migrant and the large-scale elevational migrant differ in their capacity to regulate the rate and efficiency of glycolytic enzyme activity.

Elevational migrants exhibited different levels of plasticity in lactate metabolism

We observed opposite patterns of protein level variation for two subunits of lactate dehydrogenase (LDH) in the flight muscles of the two species (Fig. 4d). LDH is the terminal enzyme of anaerobic glycolysis and consists of two main subunit types: A (previously called M for muscle) and B (previously called H for heart)^35,36. LDH-A metabolizes the conversion of pyruvate to lactate, and LDH-B catalyzes the conversion of lactate to pyruvate. These two subunits thus regulate cellular pyruvate and lactate homeostasis³⁵. We found that LDH-A exhibited greater protein expression plasticity across the four acclimation treatments in the small-scale migrant compared to the large-scale migrant. Conversely, LDH-B showed greater protein expression plasticity in the large-scale migrant than in the small-scale migrant (Fig. 4b, c). This suggests that the flight muscle of the small-scale elevational migrant has greater plasticity in lactate production but less plasticity in utilizing lactate as an aerobic substrate compared to the large-scale elevational migrant.

We further explored the dynamics of the two LDH subunits across different experimental treatments for the elevational migrants. In response to temperature change, LDH-B exhibited low protein levels in both species under warm conditions and high levels under cold conditions. However, statistical significance was observed only in the large-scale elevational migrant (t-test with FDR adjusted P < 0.05, Fig. 4d). In contrast, LDH-A showed low protein expression in cold conditions and high expression in warm conditions for both species (large-scale elevational migrant, FDR adjusted P < 0.001; small-scale elevational migrant, FDR adjusted P = 0.06, Fig. 4d). Regarding PO₂ change, no significant differences were observed in LDH-A and LDH-B expressions between hypoxic and normoxic conditions for either species (FDR adjusted P > 0.05, Fig. 4d). Collectively, these dynamics suggest that modulation of LDH-A and LDH-B subunits is primarily driven by temperature rather than hypoxia. This points to an autoregulatory loop in which activation of LDH-B leads to the reciprocal suppression of LDH-A.

Discussion

The mechanisms underlying phenotypic plasticity that co-varies with the seasonal elevational ranges of migratory species remain largely unknown, yet they are critical to understanding how organisms respond to rapid environmental change. Through common-garden experiments, multi-tissue gene expression analysis, and flight muscle proteomic analysis of two closely related species with different seasonal elevational movements, we find that phenotypic plasticity is greater in response to temperature change in the species that experiences higher seasonal temperature variation (i.e., T. indicus, the small-scale elevational migrant). In contrast, the species that experiences greater seasonal variation in PO₂ exhibited greater plasticity to changes in PO₂ (i.e., T. chysaeus, the large-scale elevational migrant). The large-scale elevational migrant, which migrates across a broad elevational range (i.e., approximately 500–4600 m a.s.l.), benefits from enhanced plasticity to cope with variable PO₂. In contrast, the small-scale elevational migrant, which migrates within a narrower elevational range (2000-4200 m a.s.l.)^22,23, requires greater plasticity to respond to seasonal temperature fluctuation. Thus, each species’ plasticity appears to be shaped by the selective pressures of their distinct seasonal elevational ranges and migratory behaviors.

Our study also provides insight into tissue-specific responses to changes in temperature and PO₂. For example, the large-scale elevational migrant exhibits greater gene expression plasticity in heart tissue compared to the small-scale elevational migrant. This difference could be attributed to the heart’s high oxidative demand^37,38, which requires greater plasticity to cope with PO₂ variation across a wide range of elevation. Additionally, we observed distinct plastic responses in flight muscle glycolytic activity between the two species. The flight muscle of bush-robins consists of two types of muscle fibers: fast oxidative fibers and fast glycolytic fibers²⁴. The fast oxidative fibers are primarily responsible for sustaining aerobic activity, while the fast glycolytic fibers support anaerobic metabolism^39,40. Our previous research has shown that while both species have similar muscle mass of fast oxidative fibers, the small-scale elevational migrant possesses a greater muscle mass of fast glycolytic fibers than the large-scale elevational migrant²⁴.

This observation may reflect a greater need for the small-scale elevational migrant to modify flight muscle glycolytic energy production. The small-scale elevational migrant remains at higher elevation throughout the year and must cope with large fluctuation in temperature compared to the large-scale elevational migrant. Glycolytic fibers contribute little to sustained thermogenesis and may even be detrimental during colder periods of the year⁴¹. Our previous work has shown that the small-scale elevational migrant has increased glycolytic fiber mass in warmer months, potentially benefiting from a reduction in glycolytic fiber mass and anaerobic activity during cold months²⁴. Alternatively, the observed differences may be attributed to the distinct migration behaviors of the two species. While glycolytic fibers can support high-intensity burst flight, their low mitochondrial content makes them ill-suited for long durations of sustained flight, which is necessary for large distance movements (i.e., large-scale migrant). Regardless of the potential causes, our findings suggest that glycolytic fibers in elevational migrants represent a unique functional context, where variation in muscle phenotypes is intricately linked to the life histories of elevational migrants.

Both elevational migrants also demonstrate temperature-driven alternations in LDH-A and LDH-B subunits. However, the plasticity of these subunits differs between the two species, potentially aiding in the removal of lactate from muscle tissue or modulating glycolytic activity in line with their distinct migratory behaviors^42,43. The large-scale elevational migrant shows greater plasticity in LDH-B, which may allow it to adjust lactate levels during both up-slope and down-slope migration. This plasticity could help the large-scale elevational migrant avoid flight muscle fatigue resulting from pH imbalance during longer migratory movements. The small-scale elevational migrant, however, which stays at higher elevations year-round and faces harsher cold conditions, exhibits greater plasticity in LDH-A within its flight muscle. This may help maintain energy balance by efficiently synthesizing ATP, ensuring a better match between ATP demand and ATP supply^44,45. Interestingly, while LDH-A and LDH-B alternation has typically been viewed as a response to hypoxia^46,47,48,49, we find that temperature-driven alternation is more pronounced than hypoxia-driven alternation in these species.

In summary, phenotypic plasticity plays a crucial role in helping animals cope with extreme climatological events⁵⁰. However, the extent to which plastic responses to changes in temperature and PO₂ co-vary with seasonal elevational distributions in montane fauna remains poorly understood. This knowledge gap is becoming increasingly critical in light of rapid global change, particularly in mountain regions, which harbor high levels of biodiversity^51,52,53,54. By combing common-garden experiments with gene and protein expression analyses, our study demonstrates that the phenotypic plasticity of the two studied elevational migrants is largely shaped by the seasonally environmental fluctuations they experience annually. The small-scale elevational migrant exhibited greater plasticity in response to temperature change, while the large-scale elevational migrant exhibited greater plasticity in response to PO₂ change. Furthermore, when examining plastic responses to changes in temperature and PO₂ independently, both species exhibited greater plasticity in response to temperature change than PO₂ change. With the ongoing impacts of global climate change, temperatures at higher elevation are expected to increase, while PO₂ is likely to remain stable. These findings enhance our understanding of how montane organisms may adapt to their changing environments.

Methods

Sampling sites

We captured the small-scale elevational migrant (T. indicus, n = 21) and the large-scale elevational migrant (T. chrysaeus, n = 23) from wild, sympatric populations in Hailuogou Valley, Sichuan, China (29.58° N, longitude, 101.99° E). All birds were captured between 3000 and 3400 m a.s.l. as they moved upslope to breeding grounds from April to early May. Only males in definitive adult plumage were included in this study, as we previously identified substantial physiological differences within bush-robin species across age and sex classes²⁴. To ensure sufficient sample sizes, experiments were conducted over two years (2015 and 2016). Prior to experiments, birds were held in captivity at 3000 m a.s.l. in individual mesh enclosures for 5–35 days. During this period, birds acclimated to captivity under ambient environment conditions before being transported to Sichuan University in Chengdu (~500 m a.s.l.), where acclimation experiments were conducted. All birds were collected with authorization from the Gongga Shan National Nature Reserve and the Hailuogou Forestry Bureau. Birds were handled in accordance with the Institutional Animal Care and Use Committee at the University of Chicago (IACUC Protocol #72394).

Acclimation experiment

We exposed bush-robins to simulated environmental conditions using a two-by-two experimental design, manipulating temperature and PO₂. Four to six individuals per species were randomly assigned to one of four experimental treatments: (1) warm temperature and normoxia, (2) cold temperature and normoxia, (3) warm temperature and hypoxia, and (4) cold temperature and hypoxia (warm temperature = 27°C ± 2°C; cold temperature = 8 ± 2 °C; normoxia = 151 mmHg; hypoxia = 92 mmHg). Warm treatments corresponded to a general thermoneutral zone in which adult birds can maintain body temperature with minimal energy expenditure, while cold treatments corresponded to temperatures below the general thermalneutral zone for birds. Both cold and warm treatments fell within the range of temperatures experienced by the two species in the wild.

Normoxia and hypoxia treatments corresponded to PO₂ levels at approximately 500 and 4650 m a.s.l., respectively. The latter corresponds to the upper limit recorded for the two species; the former represents the lower elevational limit for the large-scale elevational migrant. Although the small-scale migrant’s lower elevational limit is 2000 m a.s.l., we set 500 m a.s.l. as the lower limit for both species to ensure consistent environmental conditions across treatments and to maximize the potential differences in environmental conditions for testing species’ responses. The hypoxic group was exposed to an atmosphere simulating the partial pressure of oxygen at 4650 m a.s.l. (92 mmHg, 58% of 160 mmHg, the PO₂ at sea level). We used vacuum pumps to alter oxygen pressure for the hypoxic treatment groups. The normoxia group was acclimated to PO₂ conditions at 500 m a.s.l. (151 mmHg, 94% of 160 mmHg, the PO₂ at sea level), corresponding to the elevation of Chengdu, where the acclimation experiments were conducted. PO₂ values were calculated based on experimental atmospheric pressure × 21%, representing estimated PO₂ levels.

Photoperiod was kept constant to match the ambient photoperiod for Hailuogou Valley during the experimental period (13.5 light hours/10.5 dark hours). Birds were housed individually in experimental chambers without visual access to other individuals.

Upon arrival to Sichuan University, birds were given a two-day buffer period to recover from any stress induced by transport. After this adjustment period, birds were placed into their experimental treatments. For the first two days of the experiment, birds in hypoxia treatments were kept at 113 mmHg (equivalent of 3000 m a.s.l.) to reduce any initial shock and potential injury, after which PO₂ was reduced to 92 mmHg, the full hypoxia treatment. Birds acclimated to their full treatment environments for two weeks. The two-week acclimation period was determined based on the birds’ ecological characteristics and logistical considerations. First, bush robins at our field sites typically take a few days to move between breeding and nonbreeding elevations, indicating their ability to shift elevations rapidly²³. Second, acclimation experiments are logistically challenging for small wild-caught birds due to high mortality rates in captivity. We selected a two-week acclimation period to maximize animal welfare.

Birds in cold and hypoxia treatments were removed from their treatment environments for less than 15 min every two days for feeding and chamber cleaning. Birds in warm and normoxia treatments experienced constant conditions during feeding and chamber cleaning, as these conditions matched the ambient condition outside the experimental chambers. After the two-week experimental treatment, we collected four tissues: heart (T. chrysaeus, n = 23; T. indicus, n = 21), flight muscle (T. chrysaeus, n = 23; T. indicus, n = 21), liver (T. chrysaeus, n = 14; T. indicus, n = 14), and lung (T. chrysaeus, n = 22; T. indicus, n = 21). Tissues were all collected within 20 min of euthanasia to ensure that the birds experienced the same treatment time. The samples were immediately flash-frozen in liquid nitrogen, and transferred to a −80 °C freezer, where they remained until RNA sequencing.

RNA extraction, sequencing, expression quantification and normalization

In 2022, we conducted RNA sequencing for four tissue types: heart, flight muscle, liver, and lung (Supplementary Table S1). RNA was extracted using the RNA Easy Fast Tissue/Cell Kit DP451 (Tiangen Biotech, China), following the manufacturer’s protocol. RNA-seq libraries were constructed according to the Illumina manufacturer’s protocol. In total we sequenced 159 transcriptomes across the four tissues for the two species (Supplementary Data 1). RNA sequencing was performed using 150 bp paired-end reads using an Illumina Novaseq 6000 platform (Berry Genomics, China). After filtering low-quality reads, adapter-contaminated reads, and reads with >10% N content using FastQC v.0.11⁵⁵ and Fastp v.0.23.4⁵⁶, a total of 3,887,131,776 cleaned reads (Supplementary Table S2) were mapped to the genome of the Collared flycatcher (Ficedula albicollis, https://www.ncbi.nlm.nih.gov/datasets/genome/GCF_000247815.1/) using STAR v. 2.7.9.a⁵⁷. Gene expression levels were quantified using RSEM v. 1.3.0⁵⁸. We retained only genes expressed in 80% of the samples with TPM > 1 in every sample, resulting in 11,357 genes for T. chrysaeus and 11,447 genes for T. indicus. This yielded a dataset of 11,077 one-to-one orthologs between the two species, which was used for subsequent analyses.

Since the common garden experiments were conducted over two years (2015 and 2016), we compared the median expression levels of all orthologous genes between years for the heart, flight muscle, and lung (liver was only collected in 2016) using a Wilcoxon rank sum exact test. To assess whether the duration of captivity influenced gene expression patterns, we divided the captivity period into five time windows of seven days each and compared the median expression levels of all orthologous genes between pairs of groups using a Wilcoxon rank sum exact test with FDR correction multiple comparisons. No batch effects were detected between years (P > 0.05, Supplementary Fig. S13) or between individuals held in captivity for different lengths of time (P > 0.05, Supplementary Fig. S14). Therefore, data from both experimental years were pooled for subsequent analyses.

To normalize gene expression between the small-scale and large-scale elevational migrants, we identified 1000 conserved genes with the lowest coefficient of variance and calculated their median expression levels. These values were used to compute scaling factors, which were then applied to normalize gene expression across all samples^59,60.

Genetic variation across acclimation treatments and sampling years

To determine whether gene and protein expression differences between acclimation treatments resulted from plasticity rather than genetic variation, we quantified genetic variation (nucleotide diversity and genetic heterozygosity) among the four acclimation treatments and between the two sampling years. We used a Mann-Whitney U test with FDR correction for multiple comparisons (acclimation treatments) to test for genetic difference within treatments or years.

Identifying differentially expressed genes

We used DESeq2 v.1.22.1⁶¹ to identify differentially expressed genes (DEGs) between the following comparisons: cold/hypoxia vs. warm/normoxia (to assess responses to combined changes in temperature and PO₂), cold/hypoxia vs. warm/hypoxia and cold/normoxia vs. warm/normoxia (to assess responses to changes in temperature), and cold/hypoxia vs. cold/normoxia and warm/hypoxia vs. warm/normoxia (to assess responses to changes in PO₂), for the two species, respectively. DEGs were defined as those with |log₂FoldChange | ≥ 1 and false discovery rate (FDR) adjusted P-value < 0.05. A two-tailed exact binomial test was used to determine whether the number of DEGs differed significantly between the two species or between acclimation treatments. The code for identifying DEGs is provided in Supplementary code: R-script for DEG.

Co-expression network analysis

We used WGCNA v. 1.41-1²⁹ to identify regulatory modules correlated with acclimation treatments. Two levels of WGCNA comparisons were performed. The first comparison included all four acclimation treatments for each species, enabling the identification of hub genes suitable for comparing gene expression plasticity arising from combined changes in temperature and PO₂ versus changes in temperature and PO₂ independently. The second WGCNA comparison focused on interspecific differences in response to changes in temperature and PO₂ independently. For this analysis, we conducted two independent WGCNA analyses: one identifying hub genes associated with cold/hypoxia and warm/hypoxia conditions (changes in temperature) and another identifying hub genes associated with cold/hypoxia and cold/normoxia conditions (changes in PO₂). The approach minimized confounding effects from other acclimation conditions.

For each tissue, we calculated Pearson correlations between each pair of genes and constructed regulatory networks. A soft threshold method was used to compute an adjacency matrix, ensuring an approximately scale-free topology. Topological overlap was used to create cluster dendrograms based on hierarchical clustering. Modules were identified as branches of the cluster tree using the dynamic tree-cutting method. To identify regulatory modules associated with acclimation conditions, we performed principal component analysis (PCA) to summarize gene expression modules using the blockwiseModules function in WGCNA (minimum module size = 30, minimum height for merging modules = 0.25, with other parameters set to default). The first principal component (PC1) represented overall module expression, and module eigengene values of PC1 were used to test for associations between module expression and acclimation conditions using Pearson correlation (via the “cor” function). P values were determined using a Student’s asymptotic test (via the “corPvalueStudent” function), with significant regulatory modules identified at a P value < 0.05. Hub genes (genes strongly associated with acclimation treatments) were identified as those with P < 0.05. The code for WGCNA is provided in Supplementary code: R-script for WGCNA.

Assessment of module robustness

Due to the smaller sample size for liver tissue (n = 13, compared to the required n > 15 for WGCNA³³), we evaluated robustness of liver analysis using jackknife resampling. Specifically, 92% of the individuals (12 out of 13) were randomly selected, and the WGCNA procedure was repeated 13 times. We assessed module robustness using two approaches. First, correlation between module eigengene scores from each jackknife resample and the full dataset was evaluated by Jaccard index and Dice similarity coefficient. Second, we calculated vector analysis on the jackknife resampling modules, and compared vector lengths between the large-scale and small-scale elevational migrants. Vector analyses of the jackknife resampling modules yielded results consistent with those from the full data (Fig. S10).

Vector analysis

Based on the hub genes identified by WGCNA, we conducted vector analyses as described in Refs. ^31,32 to quantify gene expression plasticity between acclimation treatments for the two species. Specifically, we used PCA to reduce the large hub gene expression dataset into a smaller number of quantitative axes representing expression differentiation. We retained individual scores on the principal component axes that explained > 85% of the gene expression differentiation and calculated the centroid for each treatment on these axes to obtain vectors in gene expression space. These vectors were used to summarize multivariate mean values of gene expression and quantify plasticity in terms of vector lengths and angles.

Vector length was calculated as the multivariate Euclidean distance between centroids of vectors for different acclimation treatments. Using the vectors of the cold/hypoxia condition as an approximate baseline (i.e., the environmental conditions corresponding to the breeding elevations of the two elevational migrants), we calculated vector lengths between cold/hypoxia and each of the warm/hypoxia, cold/normoxia, and warm/normoxia for each species (Fig. 3a). Vector lengths thus represent the magnitude of gene expression difference in response to changes in temperature (cold/hypoxia vs. warm/hypoxia), PO₂ (cold/hypoxia vs. cold/normoxia) and a combination of temperature and PO₂ (cold/hypoxia vs. warm/normoxia).

Vector angles were calculated using the dot product of vectors, which corresponds to the arc-cosine of the Pearson correlation between vectors. Using the cold/hypoxia condition as the baseline (i.e., a standardized reference vector), we calculated the vector angle between cold/normoxia and warm/normoxia to assess gene expression difference in response to changes in temperature, and the vector angle between warm/normoxia and warm/hypoxia to assess gene expression differences in response to changes in PO₂ (Fig. 3b). Vector angles measure the direction of divergence in gene expression space and quantified the relative contribution of different gene expressions to divergence resulting from changes in temperature and PO₂ independently.

To generate effect sizes for vector angles and lengths, we used Monte Carlo simulation of P-matrices with function “rnorm” in R to generate 5000 simulated datasets based on the mean and covariance of vector angle and length, respectively^62,63. This approach accounts for sampling variance across different matrices. Using these distributions, we tested whether vector angles and lengths differed between small-scale and large-scale elevational migrants using a Wilcoxon test with FDR correction for multiple comparisons. The code for vector analyses is provided in Supplementary code: R-script for vector analysis.

Functional enrichment analysis

We performed functional enrichment analysis using DAVID (Database for Annotation, Visualization and Integrated Discovery 2021, DAVID Knowledgebase, v2023q4, https://david.ncifcrf.gov/)⁶⁴, with chicken (Gallus gallus) as the reference background. Significantly overrepresented GO terms and pathways were identified by an FDR adjusted P-value threshold of < 0.1, with 13,534 genes from chicken as the background get set (default in the DAVID database). The significantly enriched GO terms were categorized according to their biological functions in the Biological Process category.

Gene expression plasticity in metabolic pathways

We specifically evaluated gene expression plasticity for the hub genes involved in metabolic pathways, as the largest proportions of these genes were enriched for metabolic activities (Supplementary Fig. S11). We measured the coefficient of variation (CV) of gene expression across the four acclimation treatments for hub genes in fatty acid metabolism, glycolysis/gluconeogenesis, and the citrate cycle pathways. We tested whether gene expression plasticity in specific metabolic pathways differed between the two species using a Mann-Whitney U test and a significant level of P < 0.05.

Flight muscle proteomic analysis

We conducted flight muscle proteomic analysis to investigate glycolytic activities in the two species. Flight muscle proteins were extracted and digested as described in ref. ⁶⁵. Each protein sample was diluted in 200 μl phosphate-buffered saline (0.01 mol/L) containing a protease inhibitor (50×, Sigma, USA). The diluted samples were centrifuged 12,000 × g at 4 °C for 15 min, and the supernatant was retained for further analysis. Proteins in the supernatant were precipitated with 100 μl 100 mM TEAB and digested with trypsin (PROMEGA, protein-to-trypsin ratio = 1:50) overnight at 37 °C. The trypsin-digested proteins were desalted using a BCA protein Assay Kit (Pierce, Rockford, IL, USA). Purified proteins were analyzed using a data-independent acquisition (DIA) proteomics approach. DIA quantitative analyses were performed using a liquid chromatography–mass spectrometry (LC–MS) system consisting of a UPLC M-Class platform (Waters, USA) coupled with a Q Exactive HF MS spectrometer (Thermo Fisher, USA). We used an elution model for LC and the setting mode for MS as described in ref. ⁶⁶. Protein identification was refined using iRT (AB SCIEX, USA) calibration, with stringent criteria (FDR) for protein and peptide identification to ensure high confidence in the quantitative analysis.

We measured the CV of protein expression across four acclimation treatments for each identified hub gene in the glycolysis/gluconeogenesis pathway for the two species. For each CV, we used a bootstrap approach with 1000 permutations to generate 95% confidence interval^67,68. We used the Kolmogorov-Smirnov (K-S) test to assess statistical differences in protein expression CVs between the two species.

Data availability

The transcriptional data from Tarsiger indicus and Tarsiger chrysaeus generated in this study have been deposited in the Genome Sequence Archive (https://ngdc.cncb.ac.cn/gsa/) in the National Genomics Data Center under the accession number CRA024711. Other data can be found in Supplementary Data files 1–3. Source data are provided with this paper.

Code availability

Analysis scripts can be found in Github (https://github.com/xueboning/Himalayan-Bush-robins) and the supplementary codes⁶⁹.

References

Berthold P., Gwinner, E. & Sonnenschein, E. Avian Migration (Springer Science & Business Media, 2013).
Newton, I., The Migration Ecology of Birds (Elsevier Science, 2023).
Åkesson, S. & Hedenström, A. How migrants get there: migratory performance and orientation. BioScience 57, 123–133 (2007).
Article Google Scholar
Dingle, H., Migration: The Biology of Life on the Move (Oxford University Press, 2014).
Winkler, D. W. et al. Cues, strategies, and outcomes: how migrating vertebrates track environmental change. Mov. Ecol. 2, 10 (2014).
Article Google Scholar
Boyle, W. A., Sandercock, B. K. & Martin, K. Patterns and drivers of intraspecific variation in avian life history along elevational gradients: a meta-analysis. Biol. Rev. 91, 469–482 (2016).
Article Google Scholar
Hsiung, A. C., Boyle, W. A., Cooper, R. J. & Chandler, R. B. Altitudinal migration: Ecological drivers, knowledge gaps, and conservation implications. Biol. Rev. 93, 2049–2070 (2018).
Article PubMed Google Scholar
West, J. B. Prediction of barometric pressures at high altitudes with the use of model atmospheres. J. Appl. Physiol. 81, 1850–1854 (1996).
Article CAS PubMed Google Scholar
Ivy, C. M. & Williamson, J. L. On the physiology of high-altitude flight and altitudinal migration in birds. Integr. Comp. Biol. 3, icae062 (2024).
Google Scholar
Williamson, J. L. et al. Extreme elevational migration spurred cryptic speciation in giant hummingbirds. Proc. Natl Acad. Sci. USA. 121, e2313599121 (2024).
Article CAS PubMed PubMed Central Google Scholar
Meir, J. U. & Milsom, W. K. High thermal sensitivity of blood enhances oxygen delivery in the high-flying bar-headed goose. J. Exp. Biol. 216, 2172–2175 (2013).
CAS PubMed Google Scholar
Kisdi, É Dispersal: Risk spreading versus local adaptation. Am. Nat. 159, 579–596 (2002).
Article PubMed Google Scholar
Storz, J. F., Scott, G. R. & Cheviron, Z. A. Phenotypic plasticity and genetic adaptation to high-altitude hypoxia in vertebrates. J. Exp. Biol. 213, 4125–4136 (2010).
Article PubMed PubMed Central Google Scholar
Bonamour, S., Chevin, L. M., Charmantier, A. & Teplitsky, C. Phenotypic plasticity in response to climate change: the importance of cue variation. Philos. Trans. R. Soc. B 374, 20180178 (2019).
Article Google Scholar
Kelly, M. Adaptation to climate change through genetic accommodation and assimilation of plastic phenotypes. Philos. Trans. R. Soc. Lond. B Biol. Sci. 374, 20180176 (2019).
Article PubMed PubMed Central Google Scholar
Brennan, R. S. et al. Loss of transcriptional plasticity but sustained adaptive capacity after adaptation to global change conditions in a marine copepod. Nat. Commun. 13, 1147 (2022).
Article ADS CAS PubMed PubMed Central Google Scholar
Menon, T., Ramesh, V. & Barve, S. Himalayan birds that show the greatest elevational shifts remain within the narrowest thermal regimes. Glob. Ecol. Biogeogr. 32, 2111–2121 (2023).
Article Google Scholar
Anatskaya, O. V., Sidorenko, N. V., Beyer, T. V. & Vinogradov, A. E. Neonatal cardiomyocyte ploidy reveals critical windows of heart development. Int. J. Cardiol. 141, 81–91 (2010).
Article PubMed Google Scholar
Marchese, C. Biodiversity hotspots: A shortcut for a more complicated concept. Glob. Ecol. Conserv. 3, 297–309 (2015).
Google Scholar
Barve, S., Dhondt, A. A., Mathur, V. B. & Cheviron, Z. A. Life-history characteristics influence physiological strategies to cope with hypoxia in Himalayan birds. Proc. R. Soc. B. 283, 1843 (2016).
Article Google Scholar
Dixit, S., Joshi, V. & Barve, S. Bird diversityy of the Amrutganga Valley, Kedarnath, Uttarakhand, India with an emphasis on the elevational distribution of species. Check List 12, 1–11 (2016).
Article Google Scholar
Grimmett, R. et al. A Guide to the Birds of India, Pakistan, Nepal, Bangladesh, Bhutan, Sri Lanka, and the Maldives. https://www.google.com/books/edition/A_Guide_to_the_Birds_of_India_Pakistan_N/fEp4QgAACAAJ?hl=en (1999).
Collar, N. Golden Bush-Robin (Tarsiger chrysaeus), White-browed Bush-Robin (Tarsiger indicus), version 1.0. In Birds of the World (eds del Hoyo, J. et al.) (Cornell Lab of Ornithology, 2020).
DuBay, S. G. et al. Life history predicts flight muscle phenotype and function in birds. J. Anim. Ecol. 89, 1262–1276 (2020).
Article PubMed Google Scholar
Beaman, J. E., White, C. R. & Seebacher, F. Evolution of plasiticty, mechanistic link between development and reversible acclimation. Trends Ecol. Evol. 31, 237–249 (2016).
Article PubMed Google Scholar
Coolon, J. D., McManus, C. J., Stevenson, K. R., Graveley, B. R. & Wittkopp, P. J. Tempo and mode of regulatory evolution in Drosophila. Genome Res. 24, 797–808 (2014).
Article CAS PubMed PubMed Central Google Scholar
Hodgins-Davis, A., Rice, D. P. & Townsend, J. P. Gene expression evolves under a house-of-cards model of stabilizing selection. Mol. Biol. Evol. 32, 2130–2140 (2015).
Article CAS PubMed PubMed Central Google Scholar
Pollock, H. S., Brawn, J. D., Agin, T. J. & Cheviron, Z. A. Differences between temperate and tropical birds in seasonal acclimatization of thermoregulatory traits. J. Avian Biol. 50, e02067 (2019).
Article Google Scholar
Langfelder, P. & Horvath, S. WGCNA: an R package for weighted correlation network analysis. BMC Bioinform. 9, 559 (2008).
Article Google Scholar
She, H. S. et al. Gene expression plasticity followed by genetic change during colonization in a high-elevation environment. Elife 12, RP86687 (2024).
Article PubMed PubMed Central Google Scholar
Adams, D. C. & Collyer, M. L. A general framework for the analysis of phenotypic trajectories in evolutionary studies. Evolution 63, 1143–1154 (2009).
Article PubMed Google Scholar
Stuart, Y. E. et al. Contrasting effects of environment and genetics generate a continuum of parallel evolution. Nat. Ecol. Evol. 22, 158 (2017).
Article Google Scholar
Ballouz, S., Verleyen, W. & Gillis, J. Guidance for RNA-seq co-expression network construction and analysis: safety in numbers. Bioinformatics 31, 2123–2130 (2015).
Article CAS PubMed Google Scholar
Valladares, F., Sanchez-Gomez, D. & Zavala, M. A. Quantitative estimation of phenotypic plasticity: bridging the gap between the evolutionary concept and its ecological applications. J. Ecol. 94, 1103–1116 (2006).
Article Google Scholar
Storey, K. B. Comparative enzymology-new insights from studies of an “old” enzyme, lactate dehydrogenase. Comp. Biochem. Physiol. 199, 13–20 (2016).
Article ADS CAS Google Scholar
Dawson, D. M., Kaplan, N. O. & Goodfriend, T. L. Lactic dehydrogenases: Functions of the two types rates of synthesis of the two major forms can be correlated with metabolic differentiation. Science 143, 929–933 (1964).
Article ADS CAS PubMed Google Scholar
Grynberg, A. & Demaison, L. Fatty acid oxidation in the heart. J. Cardiovasc. Pharmacol. 28, 11–17 (1996).
Google Scholar
Giordano, F. J. Oxygen, oxidative stress, hypoxia, and heart failure. J. Clin. Invest. 115, 500–508 (2005).
Article CAS PubMed PubMed Central Google Scholar
Meyers, R. A. & Stakebake, E. F. Anatomy and histochemistry of spread-wing posture in birds. 3. Immunohistochemistry of flight muscles and the “shoulder lock” in albatrosses. J. Morphol. 263, 12–29 (2005).
Article PubMed Google Scholar
Dakin, R., Segre, P. S., Straw, A. D. & Altshuler, D. L. Morphology, muscle capacity, skill, and maneuvering ability in hummingbirds. Science 359, 653–657 (2018).
Article ADS CAS PubMed Google Scholar
Marsh, R. L. & Dawson, W. R. Avian adjustments to cold. In Animal Adaptation to Cold. 205–253 (Springer, 1989).
Banerjee, S. & Chaturvedi, C. M. Migratory preparation associated alterations in pectoralis muscle biochemistry and proteome in Palearctic-Indian emberizid migratory finch, red-headed bunting, Emberiza bruniceps. Comp. Biochem Physiol. Part D. Genomics Proteom. 17, 9–25 (2016).
Article CAS Google Scholar
Sheafor, B. A. Metabolic enzyme activities across an altitudinal gradient: an examination of pikas (genus Ochotona). J. Exp. Biol. 206, 1241–1249 (2003).
Article CAS PubMed Google Scholar
McFarlan, J. T., Bonen, A. & Guglielmo, C. G. Seasonal upregulation of fatty acid transporters in flight muscles of migratory white-throated sparrows (Zonotrichia albicollis). J. Exp. Biol. 212, 2934–2940 (2009).
Article CAS PubMed Google Scholar
Zajac, D. M., Cerasale, D. J., Landman, S. & Guglielmo, C. G. Behavioral and physiological effects of photoperiod-induced migratory state and leptin on Zonotrichia albicollis: II. Effects on fatty acid metabolism. Gen. Comp. Endocrinol. 174, 269–275 (2011).
Article CAS PubMed Google Scholar
Dezwaan, A. & Wijsman, T. C. M. Anaerobic metabolism in Bivalvia (Mollusca) - characteristics of anaerobic metabolism. Comp. Biochem. Physiol. B 54, 313–323 (1976).
Article CAS Google Scholar
Collicutt, J. M. & Hochachka, P. W. The anaerobic oyster heart: coupling of glucose and aspartate fermentation. J. Comp. Physiol. 115, 147–157 (1977).
Article CAS Google Scholar
Saz, H. J. Energy metabolisms of parasitic helminths: adaptations to parasitism. Ann. Rev. Physiol. 43, 323–341 (1981).
Article CAS Google Scholar
Hochachka, P. W., Stanley, C., McKenzie, D. C., Villena, A. & Monge, C. Enzyme mechanisms for pyruvate-to-lactate flux attenuation: study of Sherpas, Quechuas, and hummingbirds. Int. J. Sports Med. 13, S119–S122 (1992).
Article CAS PubMed Google Scholar
Regan, C. E. & Sheldon, B. Phenotypic plasticity increases exposure to extreme climatic events that reduce individual fitness. Glob. Chang. Biol. 29, 2968–2980 (2023).
Article CAS PubMed PubMed Central Google Scholar
Fjeldså, J., Bowie, R. C. K. & Rahbek, C. The role of mountain ranges in the diversification of birds. Annu. Rev. Ecol. Evol. Syst. 43, 323–341 (1981).
Google Scholar
Freeman, B. G., Scholer, M. N., Ruiz-Gutierrez, V. & Fitzpatrick, J. W. Climate change causes upslope shifts and mountaintop extirpations in a tropical bird community. Proc. Natl Acad. Sci. USA 115, 11982–11987 (2018).
Article ADS CAS PubMed PubMed Central Google Scholar
Inouye, D. W., Barr, B., Armitage, K. B. & Inouye, B. D. Climate change is affecting altitudinal migrants and hibernating species. Proc. Natl Acad. Sci. USA 97, 1630–1633 (2000).
Article ADS CAS PubMed PubMed Central Google Scholar
Chen, Y. L. et al. The combination of genomic offset and niche modelling provides insights into climate change-driven vulnerability. Nat. Commun. 13, 4821 (2022).
Article ADS CAS PubMed PubMed Central Google Scholar
Andrews, S. FastQC: A Quality Control tool for High Throughput Sequence Data. http://www.bioinformatics.babraham.ac.uk/projects/fastqc/ (2014).
Chen, S. F., Zhou, Y. Q., Chen, Y. R. & Gu, J. fastp: an ultra-fast all-in-one FASTQ preprocessor. Bioinformatics 34, 884–890 (2018).
Article Google Scholar
Dobin, A. et al. STAR: ultrafast universal RNA-seq aligner. Bioinformatics 29, 15–21 (2013).
Article CAS PubMed Google Scholar
Li, B. & Dewey, C. N. RSEM: accurate transcript quantification from RNA-Seq data with or without a reference genome. BMC Bioinform. 12, 323 (2011).
Article CAS Google Scholar
Hao, Y. et al. Comparative transcriptomics of 3 high-altitude passerine birds and their low-altitude relatives. Proc. Natl. Acad. Sci. USA 116, 11851–11856 (2019).
Article ADS CAS PubMed PubMed Central Google Scholar
Brawand, D. et al. The evolution of gene expression levels in mammalian organs. Nature 478, 343–348 (2011).
Article ADS CAS PubMed Google Scholar
Anders, S. & Huber, W. Differential expression analysis for sequence count data. Genome Biol. 11, R106 (2010).
Article CAS PubMed PubMed Central Google Scholar
Noble, D. W. A., Radersma, R. & Uller, T. Plastic responses to novel environments are biased towards phenotype dimensions with high additive genetic variation. Proc. Natl. Acad. Sci. USA 116, 13452–13461 (2019).
Article ADS CAS PubMed PubMed Central Google Scholar
Radersma, R., Noble, D. W. A. & Uller, T. Plasticity leaves a phenotypic signature during local adaptation. Evol. Lett. 4, 360–370 (2020).
Article PubMed PubMed Central Google Scholar
Huang, D. W., Sherman, B. T. & Lempicki, R. A. Systematic and integrative analysis of large gene lists using DAVID bioinformatics resources. Nat. Protoc. 4, 44–57 (2009).
Article CAS PubMed Google Scholar
Wang, N. M. et al. Protein regulation mechanism of cold tolerance in Haemaphysalis longicornis. Insect Sci. 30, 725–740 (2023).
Article CAS PubMed Google Scholar
Agwunobi, D. O. et al. Phosphoproteomic analysis of Haemaphysalis longicornis saliva reveals the influential contributions of phosphoproteins to blood-feeding success. Front. Cell. Infect. Microbiol. 11, 769026 (2022).
Article PubMed PubMed Central Google Scholar
Hao, Y. et al. Divergent contributions of coding and noncoding sequences to initial high-altitude adaptation in passerine birds endemic to the Qinghai-Tibet Plateau. Mol. Ecol. 32, 3524–3540 (2023).
Article CAS PubMed Google Scholar
Wright, A. E. et al. Variation in promiscuity and sexual selection drives avian rate of Faster-Z evolution. Mol. Ecol. 24, 1218–1235 (2015).
Article PubMed PubMed Central Google Scholar
Xue, B. et al. Code and data for ‘Phenotypic plasticity co-varies with the elevation range of in two species of avian elevation migrants in Himalayas’. GitHub https://github.com/xueboning/Himalayan-Bush-robins (2025).

Download references

Acknowledgements

The authors thank Xueying Zhang for discussion regarding the interpretation of results. This research was funded by the National Natural Science Foundation of China (NSFCU23A20162 and NSFC32020103005 for Y.Q., and NSFC31971413 for D.L.), Third Xinjiang Scientific Expedition and Research (XJKK) (2022xjkk0205 for Y.Q.), American Museum of Natural History, American Ornithological Society, American Philosophical Society, Animal Behavior Society, National Geographic Society, Paulson Institute, and the Society for the Study of Evolution.

Author information

These authors contributed equally: Boning Xue, Huishang She, Lirong Zuo.

Authors and Affiliations

Key Laboratory of Animal Biodiversity Conservation and Integrated Pest Management, Institute of Zoology, Chinese Academy of Sciences, Beijing, China
Boning Xue, Huishang She & Yanhua Qu
College of Life Sciences, University of Chinese Academy of Sciences, Beijing, China
Boning Xue & Yanhua Qu
Hebei Key Laboratory of Animal Physiology, Biochemistry, and Molecular Biology, Hebei Collaborative Innovation Center for Eco-Environment, Ministry of Education Key Laboratory of Molecular and Cellular Biology, College of Life Sciences, Hebei Normal University, Shijiazhuang, China
Lirong Zuo & Dongming Li
Negaunee Integrative Research Center, Field Museum of Natural History, Chicago, IL, USA
John M. Bates & Shane G. DuBay
Division of Biological Sciences, University of Montana, Missoula, MT, USA
Zachary A. Cheviron
College of Life Sciences, Sichuan University, Chengdu, China
Yongjie Wu
Department of Biology, University of Texas at Arlington, Arlington, TX, USA
Shane G. DuBay
Committee on Evolutionary Biology, University of Chicago, Chicago, IL, USA
Shane G. DuBay

Authors

Boning Xue
View author publications
Search author on:PubMed Google Scholar
Huishang She
View author publications
Search author on:PubMed Google Scholar
Lirong Zuo
View author publications
Search author on:PubMed Google Scholar
John M. Bates
View author publications
Search author on:PubMed Google Scholar
Zachary A. Cheviron
View author publications
Search author on:PubMed Google Scholar
Yongjie Wu
View author publications
Search author on:PubMed Google Scholar
Dongming Li
View author publications
Search author on:PubMed Google Scholar
Shane G. DuBay
View author publications
Search author on:PubMed Google Scholar
Yanhua Qu
View author publications
Search author on:PubMed Google Scholar

Contributions

This study was conceived by S.D., Y.Q., Y.W., Z.A.C. and J.B. Common garden experiments were performed by S.D. and Y.W. The analyses were carried out by B.X., H.S., and L.Z. L.Z. and D.L. coordinated data provision. The manuscript was written by Y.Q. with contributions from S.D. and D.L.

Corresponding authors

Correspondence to Yongjie Wu, Dongming Li, Shane G. DuBay or Yanhua Qu.

Ethics declarations

Competing interests

The authors declare no competing interests.

Peer review

Peer review information

Nature Communications thanks the anonymous reviewer(s) for their contribution to the peer review of this work. A peer review file is available.

Additional information

Publisher’s note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Supplementary information

Supplementary Information

Transparent Peer Review file

Description of Additional Supplementary Files

Supplementary Data 1

Supplementary Data 2

Supplementary Data 3

Supplementary Software 1

Source data

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution-NonCommercial-NoDerivatives 4.0 International License, which permits any non-commercial use, sharing, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if you modified the licensed material. You do not have permission under this licence to share adapted material derived from this article or parts of it. The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by-nc-nd/4.0/.

Reprints and permissions

About this article

Cite this article

Xue, B., She, H., Zuo, L. et al. Phenotypic plasticity co-varies with elevational range in two avian species of elevational migrants in the Himalayas. Nat Commun 16, 5316 (2025). https://doi.org/10.1038/s41467-025-60770-w

Download citation

Received: 29 June 2024
Accepted: 30 May 2025
Published: 18 June 2025
DOI: https://doi.org/10.1038/s41467-025-60770-w