Comparison between short-term stress and long-term adaptive responses reveal common paths to molecular adaptation.

The phenotypic plasticity in responses to short-term stress can provide clues for understanding the adaptive fixation mechanism of genetic variation during long-term exposure to extreme environments. However, few studies have compared short-term stress responses with long-term evolutionary patterns; in particular, no interactions between the two processes have been evaluated in high-altitude environment. We performed RNA sequencing in embryo fibroblasts derived from great tits and mice to explore transcriptional responses after exposure to simulated high-altitude environmental stresses. Transcriptional changes of genes associated with metabolic pathways were identified in both bird and mice cells after short-term stress responses. Genomic comparisons among long-term highland tits and mammals and their lowland relatives revealed similar pathways (e.g., metabolic pathways) with that initiated under short-term stress transcriptional responses in vitro. These findings highlight the indicative roles of short-term stress in the long-term adaptation, and adopt common paths to molecular adaptation in mouse and bird cells.

Introduction

The duration of exposure to extreme environmental stress may often lead to different responses of organisms, from phenotypic to genetic adaptation, from expression shifts to adaptive mutation (Ghalambor et al., 2015; Qu et al., 2021; Qu et al., 2020). As a result, the response to short-term stress appears to be passive and plastic, while the response to long-term selection is likely fixed and irreversible. Despite separate attempts of studying short-term and long-term responses, the relationships between the two remain a debate for whether short-term plastic response may limit or promote the long-term adaptive evolution. Theories suggest that the adaptive responding ways are determined by the degree of the plasticity of short-term response: adaptive or non-adaptive (Price et al., 2003; Tang et al., 2017). If the response is adaptive, the plasticity will promote the long-term adaptive evolution. The degree of plasticity depends on how close the environmentally induced phenotype is to the local optima for the environment (Ghalambor et al., 2007; Lopez-Maury et al., 2008). For example, lizards captured in the simulated wild habitats have similar coloration with wild populations, suggesting the adaptive response to short-term stress is in the same direction as evolutionary changes (Corl et al., 2018). In contrast, non-adaptive transcriptional changes were observed in the brains of fishes raised without predators, which is in the opposite direction of adaptive changes evolved by wild fishes (Ghalambor et al., 2015). Here, we test this by comparing the long-term adaptation to extreme environments on the Qinghai-Tibet Plateau (QTP) with the expression responses to short-term stresses in two homeotherms.

The QTP is the highest plateau on the earth, and is characterized by multiple extreme environmental conditions including low temperatures, low availability of oxygen, and strong UV radiation (UVR). The high-altitude environment provides a unique opportunity to study the genetic basis of long-term adaptations. Evaluation of short-term pre-adaptation stress response often involves observations of cognitive neuroscience, morphology, and ecology (Ardila, 2016; Shah et al., 2013). Previous studies on acclimation experiment of lowland tree sparrow (Passer montanus) and rufous-collared sparrow (Zonotrichia capensis) revealed no changes in muscular phenotypes and no differences in genes expressed in muscles (Cheviron et al., 2008; Qu et al., 2020). Birds and mice at higher elevations experience colder temperatures and greater UVR exposure, although these effects are likely to be tissue-specific. However, it is not feasible to accurately compare expression patterns from different tissues of whole organisms that are subjected to these conditions. For example, UVR exposure is likely to be most important for skin, but may be less important to other tissues. Likewise, low temperatures experienced by the entire organism may affect many tissues through changes in metabolism, even if the tissues themselves remain at a constant body temperature. Indeed, many conditions are more difficult to control in short-term response experiments of whole animals subjected to different conditions compared to the interpretation of expressional differences in the isolated cell lines. Because of the availability of primary fibroblast from great tit and mice (Chen et al., 2018), we chose them as the target of artificial environment for this study about the plasticity changes.

Birds are good genomic models which can be used to investigate selection pressures that are shared with mammals. Previous studies reveal that birds and mammals have many similar physiological characteristics including isothermal body temperatures, and they are also parasitized by many common pathogens (Garant et al., 2005; Lai et al., 2013; Pap et al., 2007; Shultz and Sackton, 2019; Zhu et al., 2018). High-altitude environments can alter the oxygen-carrying capacity and energy metabolism of birds, rendering them difficult to maintain normal physiological activities (Cheviron et al., 2014). In addition, mice (Mus musculus) are important models for investigating the molecular biology and genomics of mammalian evolutionary adaptation, disease occurrence, and reproductive physiology (Chalfin et al., 2014; Emes et al., 2003; Perlman, 2016). Predictive evolutionary models of birds and mice have practical importance, for example, organisms with certain diseases respond plastically to treatments and also evolve, thus gene expression profiles used to predict plastic response to stress treatment can be utilized to estimate disease progression (Merlo et al., 2006). Genomic and transcriptomic investigations of birds and mice have identified many genes that are associated with high-altitude adaptations, but few comparative analyses have been conducted with species-specific primary cultured cells.

Here, we compare plastic expression and adaptive responses under short-term stress treatments and long-term adaptive evolution. We hypothesized that short-term stress from UVR exposure, low temperatures, and hypoxia, simulating the high-altitude extreme environments, can induce plastic gene expression responses in cells that facilitated the long-term adaptive evolution to high-altitude environments. Great tit and mouse embryonic fibroblasts were used as representative subjects for short-term stimulation experiments. Transcriptomic analyses were conducted on Parus major and M. musculus cells that were exposed to artificial high- and low-altitude environments (cold treatment, UVR exposure, and hypoxia), which was combined with whole genome resequencing data for high- vs. low-altitude P. major populations to test above hypotheses.

Results

Plastic expressional responses after short-term stimuli

We generated a total of 1.1 billion reads. The raw reads were trimmed and filtered based on their quality, and the clean reads were normalized and used for subsequent analyses (Table S1). We analyzed transcriptional patterns of great tit embryonic fibroblasts (GEF) and mouse embryonic fibroblasts (MEF) in response to short-term stimuli: cold exposure (4°C) over 1 and 4 h, hypoxia stress (3% O2) for 4,24, 0, and 1 h after UVR exposure (100 J/m2, Figures 1A and S1), and found that a total of 26 gene modules were identified in GEF (Figure 1B, magenta, green, dark gray, gray60, royalblue, red, light yellow, brown, purple, yellow, pink, light green, darkgreen, turquoise, greenyellow, blue, cyan, lightcyan, black, midnightblue, darkturquoise, salmon, dark red, tan, orange, and gray modules) and 21 gene modules were observed in MEF (Figures 1B and S2–S5 and Table S6, lightyellow, pink, gray60, turquoise, lightcyan, red, cyan, blue, tan, yellow, brown, royalblue, midnightblue, greenyellow, black, purple, magenta, lightgreen, green, salmon, and gray modules). In GEF, the most impactful environmental stimuli for gene co-expression in birds were hypoxia (blue module, cor = 0.69, p = 7e-04), followed by UVR exposure (lightcyan module, cor = −0.61, p = 0.004), with low temperatures yielding the weakest response (gray module, cor = 0.45, p = 0.05) (Figure 1B). The genes in the top five modules with highest correlation coefficients (blue, brown, turquoise, lightcyan, and tan modules) were significantly enriched in the metabolic pathways (ko01100; blue and brown modules), MAPK signaling pathway (ko04010; lightcyan module), ribosomal (ko03010; turquoise module), and spliceosome (ko03040; tan module). In MEF, genes in the top five modules with highest correlation coefficients were significantly enriched in metabolic pathways (two modules), cancer pathways, purine metabolism, and lysosome pathways. The most significant signaling pathway initiated after high-altitude external stimuli for great tit and mouse cells was the metabolic signaling pathway (Figure 1C), in which a total of 11 genes (A4GALT, ALDOC, AMPD3, ENO2, GBE1, GCLM, GMDS, LDHA, MPI, PFKL, and PGAM1) were overlapped between the great tit and mouse cells. These results suggest that metabolic pathways are the most associated with altered gene expression after high-altitude stimuli exposure.

Figure 1

Transcriptional patterns detected in GEF and MEF cell lines after high-altitude simulated conditions were significantly associated with metabolic pathways

(A) Schematic for cold, UVR exposure, and hypoxia treatment designs.

(B) The left panel shows the weighted correlation network analysis of GEFs after stimulation. Gene expression modules are shown that exhibited correlations with cold exposure over 1 (GC1) or 4 h (GC4), 4 (GH4) or 24 h (GH24) of hypoxia stress, and 0 (GU0) or 1 h (GU1) after UVR exposure. The blue and light cyan modules are specifically listed and discussed in the text. Positive correlations are indicated in red coloring, while negative correlations are indicated in blue, based on the scale to the right. The right panel shows the weighted correlation network analysis of MEFs after treatment with cold exposure over 1 (MC1) or 4 h (MC4), 4 (MH4) or 24 h (MH24) of hypoxia stress, and 0 (MU0) or 1 h (MU1) after UVR exposure. The pink and turquoise modules are highlighted and discussed in the text.

(C) The most significant pathways enriched for four modules, including two from the GEF and two from the MEF networks.

(D) The most highly abundant transcriptional factor binding motifs present within genes within the four modules of panel c (blue and light cyan modules are from the GEF profile, while the pink and turquoise modules are from the MEF profile). See also Figures S1–S5, Tables S1, S2, and S6.

The CLOVER (Cis-eLement OVERrepresentation) tool was used to identify protein motifs that were statistically overrepresented in the 200 bp upstream of each module sequence group that shared a common function. Several transcription factors were significantly enriched in these binding motifs including the ZNF384, TBP, SP1, and ETV6 motifs, blue and light cyan modules are from the GEF profile, while the pink and turquoise modules are from the MEF profile (Figure 1D and Table S7). These analyses demonstrate that responses in cells to short-term stress are likely dominated by transcriptional variation, thereby providing a genome-wide view of common and specific stress responses to hypoxia, UVR exposure, and cold temperatures in great tits and mice.

Gene set enrichment analysis (GSEA) indicated that genes significantly enriched in GEF after stress treatments were associated with oxidative phosphorylation, glycosaminoglycan degradation, arachidonic acid metabolism, and ribosomal functioning (Figures 2A and 2B). In contrast, MEF gene terms significantly enriched after stress treatments were involved in glycosaminoglycan degradation, galactose metabolism, glycolysis, starch and sucrose metabolism, lysosome functioning, nicotinate and nicotinamide metabolism, fructose and mannose metabolism, central carbon metabolism in cancer, retinol metabolism, and steroid hormone biosynthesis (Figures 2C and 2D). These results were consistent with the WGCNA analysis of GEF and MEF regulatory networks after short-term stress. These results also represent the report of stress-related activation of metabolic pathways in GEF and MEF. Further, our results provide evidence for multiple metabolic-related pathways that were significantly enriched in GEF and MEF transcriptional responses after UVR exposure, cold temperatures treatment, and hypoxia simulation.

Figure 2

Differentially expressed genes in GEF and MEF cells after short-term stimuli

(A) Gene set enrichment analysis (GSEA) plots showing representatives of the two most significantly enriched metabolic pathways in GEFs.

(B) Enrichment of KEGG signaling pathways of different gene sets from GSEA in GEFs.

(C) GSEA plots showing two representative significantly enriched pathways related to metabolism in MEFs.

(D) The ten most significantly enriched pathways comprising different gene sets from GSEA in MEFs.

(E) Abundances of differentially expressed genes (DEGs) in GEFs under different treatment with cold exposure over 1 (GC1) or 4 h (GC4), 4 (GH4) or 24 h (GH24) of hypoxia stress, and 0 (GU0) or 1 h (GU1) after UVR exposure. GC1, GC4, and GH4 groups do not have differential expression genes.

(F) Abundances of DEGs in MEFs under different treatment with cold exposure over 1 (MC1) or 4 h (MC4), 4 (MH4) or 24 h (MH24) of hypoxia stress, and 0 (MU0) or 1 h (MU1) after UVR exposure. MC1 group does not have differential expression genes.

(G) Comparison of DEG numbers between GEFs and MEFs under different treatment conditions. ∗∗∗ indicates p <0.001. Data are represented as mean ± SEM

(H) The most significantly enriched pathways of DEGs between GEFs and MEFs under different treatment conditions. See also Figures S6–S9 and Table S7.

Gene expression differences across treatments

Differentially expressed genes (DEGs) across treatments were identified in response to cold treatment, hypoxia, and UVR exposure for GEF (Figures S6 and S7). The expression profiles of DEGs after 24 h of hypoxia and 0 and 1 h after UVR exposure exhibited similar patterns in GEFs and MEFs (Figures 2E and 2F). However, DEG abundances were significantly different between GEFs and MEFs after 4 h of cold treatment, 24 h hypoxia, and 0 or 1 h after UVR exposure (Figure 2G).

The functional enrichment for DEGs showed that eight pathways were significantly enriched by GEF DEGs (corrected p <0.05) in the 24-h hypoxia-treated cells. The most significantly enriched pathway was the ECM-receptor interactions [the extracellular matrix (ECM) consists of a complex mixture of structural and functional macromolecules and plays an important role in tissue and organ morphogenesis and in the maintenance of cell and tissue structure and function] (Figure 2H and Table S7). The MAPK signaling pathway was the most significantly enriched pathway in GEF after 1 h of UVR exposure. According to our multiple sequence alignment, most of the structural features of MAPK1 (one important gene in the MAPK pathway) were found to be conserved among its homologs, such as protein kinase ATP binding site, serine/threonine-protein kinase active site, and mitogen-activated protein kinase in the S_TKc sequence. The comparison also showed that most residues in the active site were conserved in MAPK1 across the MAPK1 family (Figure S8). Functional experiments on one of the genes in the MAPK signaling pathway, MAPK1, showed that it can activate a range of cellular responses including cellular proliferation, differentiation, and apoptosis (Figure S9).

In contrast, 1 h after UVR exposure in MEF led to statistically significant pathways associated with systemic lupus erythematous, alcoholism, and viral carcinogenesis. Likewise, the most significantly enriched pathways in GEFs and MEFs were strikingly different under other treatments. For example, the HIF signaling pathway [Hypoxia-inducible factor 1 (HIF-1) is a transcription factor that functions as a master regulator of oxygen homeostasis] was the most significantly enriched signaling pathway for MEFs after 4 h of hypoxia treatment (Figure 2H). This indicates that the stress response mechanisms between mice and birds are different after short-term stimulation, and different signaling pathways are activated in response to external stimuli, suggesting a plastic response to stress in birds and mammals.

Candidate genes underlying long-term adaptations to high-altitude environments

To identify genes under positive selection (PSGs) imposed by high-altitude environment, genome-wide selective sweep analyses were conducted using XP-CLR (the cross-population composite likelihood ratio test, which is designed to run on hdf5 files representing genetic data as generated using scikit-allel). In addition, cross-population-extended haplotype homozygosity (XP-EHH) analyses were conducted based on whole genome sequencing data of 12 high-altitude and 26 low-altitude great tits (Table S2). A total of 3,277 windows were observed in the top 5% of XP-CLR scores (XPCLR score were calculated window by window across the genome, and depends on the numbers of highland and lowland population SNPs), which harbors 535 genes that were used for subsequent analyses (Figure 3A and Table S8). Four candidate genes exhibited the highest XP-CLR scores, including STEAP2, CFAP69, ZNF804B, and STEAP1 (XP-CLR score = 33.51).

Figure 3

Identification of selective sweeps in great tit genomes

(A) Genomic landscape of the XP-CLR values in the genome of highland and lowland great tit.

(B) Genomic landscape of the XP-EHH values in the genome of highland and lowland great tit.

(C) The EHHs at varying distances around the core haplotype at PARD3 and LOC107202990. See also Table S8.

XP-EHH analysis identified 2,734 candidate SNPs in highland great tits with analytical thresholds of 1.11 (66 genes, within the top 5%) and −1.1 (60 genes, within the lowest 5%), respectively (Figure 3A and Table S8). In a comparison between high- and low-altitude great tits, the top two candidate genes are PARD3 (XP-EHH value = 6.93) (Figure 3C) and LOC107202990 (XP-EHH value = 8.18) (Figure 3D and Table S8).

Relationship between stress responses and adaptive differences

Four genes were shared by the sets of DEGs and PSGs. The DEGs of GH24 included COL19A1 and SEMA5A that were also positively selected, whereas the DEGs of GU0 and GU1 showed ABCA13 and SOX5 to be overlapped with PSGs, respectively (Table S3).

In order to test whether PSGs share similar expressional profiles with DEGs, we compared (1) total connectivity (TC) of WGCNA networks, (2) level of expressional changes (i.e., fold change), and (3) intramodular characters of WGCNA networks of PSGs and DEGs. Briefly, TC measures the overall relatedness of a gene to other genes in a WGCNA network. For intramodular characters, we calculated module membership (MM) and gene significance (GS) for each gene of each module. MMs measure the correlation between the expression profile of a gene and module eigengenes, which can be considered as a representative of the gene expression profile of a module. Genes with high MM values indicate the genes are highly connected to within a module, or have high intramodular connectivity. GSs measure the correlation between a treatment and gene expression profiles. We found that PSGs always have higher TC than non-PSGs (Figure 4B), whereas the patterns in DEGs were different among treatments (Figure 4C). Within modules (for a detailed description for functional enrichment, see Table S6), we did not find consistently higher MMs in PSGs than that in non-PSGs (Figures 4F and 4G), indicating that the high TC of PSGs may be attributed to high intermodular connectivity but not high intramodular connectivity. We also found that the GSs of DEGs were almost always significantly higher than those of PSGs with a moderate to large effect size (Table S4, Figures 4F and 4G). Finally, we found different levels and directions of expressional changes in PSGs across treatments (Figures 4C and 4D). These results appear to explain the reason why the PSGs and DEGs are rarely found to be overlapped.

Figure 4

Relationships of gene expression from different cell treatments and selective sweep genes from high-altitude great tit and low-altitude great tit

(A) The overlap of the selected sweep genes in the XP-CLR and XP-EHH analysis.

(B) Connectivity of the selected genes and non-selected genes by XP-CLR and XP-EHH analysis.

(C) Connectivity between differentially expressed genes and genes with no differential expression.

(D) The expression of the selected genes and non-selected genes by XP-CLR analysis.

(E) The expression of the selected genes and non-selected genes by XP-EHH analysis.

(F) Module membership (MM) and gene significance (GS) for each gene of the blue module from the GH24.

(G) Module membership (MM) and gene significance (GS) for each gene of the brown module from the GU0. Data are represented as mean ± SEM Statistically significant associations are those marked with asterisks. ∗∗∗p <0.001, ∗∗p <0.01, ∗p <0.05. See also Figure S10, Tables S3–S5.

To evaluate whether the functional enrichment of genes within selected genomic regions of great tits was similar to the pattern of those genes observed as upregulated after short-term stress treatment of GEFs, pathway enrichment tests were conducted using the chicken genome as a reference organism. The most significantly enriched pathway comprised metabolic pathways, which was also the most significantly enriched pathway among co-expressed genes responding to short-term stress treatment (Blue Module in GEF, Pink Module in MEF, Figure 1C).

Similar functions of PSGs in high- and low-altitude mammal/bird comparisons

We tested whether the adaptive divergence observed in PSGs in the mammals also occurs in biological processes represented by Gene Ontology terms and phenotypic abnormality. We utilized g: profiler annotations to analyze the functional enrichments of PSGs in high-altitude and low-altitude mammals. We found that PSGs are generally related to cellular component, localization, metabolic process, and phenotypic abnormality (Figure 5A). According to the functional enrichment of positive selection genes and convergent genes (the independent derivation of the same gene expression pattern) in high-altitude passerines (Hao et al., 2019), especially the enrichment of KEGG pathway, we found that many functional pathways related to energy metabolism is concentrated in lipid metabolism. In order to further explore whether high-altitude passerines actually have a tendency to enrich energy metabolism, we reclassified the selected positive selection genes and convergent genes into different metabolic pathways based on the KEGG_B_class annotation, and then used the hypergeometric test to analyze these metabolic pathways. The pathway was re-enriched and analyzed, and it was found that whether it is positively selected genes or convergent genes, they are more often classified into lipid metabolism rather than carbohydrate metabolism. Specifically, for positively selected genes, 11 genes belong to lipid metabolism and seven genes belong to carbohydrate metabolism; for convergent genes, 14 genes belong to lipid metabolism and six genes belong to carbohydrate metabolism (Figure 5B). Overall, when passerines face extreme high-altitude environments, positive selection genes and convergent genes may have a significant tendency to lipid metabolism in terms of functional enrichment of energy metabolism pathways.

Figure 5

Functional enrichment analyses of PSGs in cattle, goat, pig, and six tit species

(A) GO terms are grouped into five major categories including cellular process (dark blue), cellular component (green), localization (pink), metabolic processes (purple), and phenotypic abnormality (yellow).

(B) Functional enrichments of the KEGG pathways for positive selected genes (PSGs) and convergent genes in six tit species from high- and low-altitude. The number represents the number of genes belonging to metabolism in KEGG_B_class.

Discussion

Our study used three short-term stimuli experiments to approximate high-altitude conditions (low temperature, hypoxia, and strong UVR exposure) in in vitro avian and mouse primary cells. We then compared short-term transcriptional stress responses to long-term genomic adaptive responses to investigate whether phenotypic plasticity facilitates or constrains adaptive evolution. The utilization of fibroblasts to assess responses to hypoxia, UVR exposure, or low temperature may be more reasonable. Nevertheless, these experiments were conducted because of several important considerations. There are some limitations in performing stress response studies in vitro. In contrast, most previous studies of high-elevation adaptation have focused on multi-omics data or otherwise have used morphological and physiological methods alone (Hao et al., 2019; Xiong et al., 2020). Thus, transcriptional analyses and in vitro cellular biological observations were combined in this study to evaluate the stress response mechanisms of birds and mice.

Indeed, an important focus of evolutionary biology has been the investigation of relationships between response induced by short-term ecological pressures and responses resultant from genetic evolution due to long-term transmission across generations. For example, rapid local adaptation is linked with phenotypic plasticity (Sun et al., 2020). Our results show that overlapping plastic expressional responses to short-term stress are enriched in metabolic pathways common to great tit and mouse embryonic fibroblast cells. Adaptive metabolic responses to environmental cues are intrinsic property of all cells (Davies et al., 2019). For example, exposure to cold (Stocks et al., 2004) and altitude (Bailey and Davies, 1997) increases the concentrations of circulating catecholamines. This hormonal response has metabolic implications and is notably associated with glycolysis (Gibson et al., 2017). Our results suggest that phenotypic plasticity facilitates adaptive evolution in metabolic pathways, which is consistent with the rapid metabolic adaptation to high elevation (Cheng et al., 2021; Qu et al., 2020), and with that the plasticity for color in lizards is generally adaptive and goes in the same direction as evolved changes (Corl et al., 2018). One issue that arises is that the scans for selection detect the actual targets of selection or closely linked mutations, while the transcriptional changes may be due to changes in cis (i.e. target of selection) or to changes in trans (i.e. downstream effects of changes at upstream trans-acting regulators). Future studies are needed to look at allele-specific expression patterns in heterozygous individuals as a way to identify cis-eQTL.

Previous studies of the high-altitude population exhibited an accelerated evolution of hypoxia responses (i.e., the MAPK signaling pathway) compared to lowland populations (Qu et al., 2015). Our results indicate here that birds and mammals have similar hypoxia responses after short-term stimuli and the genomic adaptive evolutionary analyses via the MAPK signaling pathway. The MAPK signaling pathway can trigger responses of cells to environmental stresses (Verma et al., 2018) and plays an important role in tumor growth and transformation, depending on angiogenesis and alterations in glucose metabolism (Folkman, 1995). Functional experiments on one of the genes, MAPK1, show that it can activate a range of cellular responses including cellular proliferation, differentiation, and apoptosis, in addition to activating immune defenses (Osborne et al., 2012; Perera et al., 2016; Santamaria and Nebreda, 2010; Upadhya et al., 2013). Future studies would be necessary to compare evolved differences in expression between high- and low-elevation birds and mammals when raised in a common environment, and ask whether they go in the same direction (the MAPK signaling pathway) or whether they generally oppose each other.

In conclusion, it is critical to further investigate short-term stress responses and longer-term adaptive evolution to high-altitude conditions in transcriptional and genomic responses of multiple species. In this study, we have compared transcriptomic profiles of cell lines from different species after short-term stress treatments. Moreover, we have demonstrated that short-term stress variations and long-term evolutionary adaptations share the common metabolic signaling pathways (Figure 6). Despite several limitations discussed above, our study suggests that phenotypic plasticity (e.g., changes in gene expression) can promote adaptive evolution, which provides valuable targets for further studies of gene expression plasticity and the complex phenotypic adaptations associated with life in extreme environments.

Figure 6

Working model for the functional enrichment of short-term stress and long-term adaptation analysis

Limitations of the study

There are two limitations to the present study. First, it is unfortunate that we do not have comparable data on high- vs. low-altitude mice to do the scans of selection. Second, our study is limited to studying the scans for selection detect the actual targets of selection or closely linked mutations, while the transcriptional changes may be due to changes in cis or to changes in trans. If transcriptional changes are in cis, we might expect to see some overlap between the stress response and the adaptive response. If transcriptional changes are in trans, there is no expectation that the adaptive targets of selection could be the same as those showing changes in expression.

STAR★Methods

Key resources table

Resource availability

Lead contact

Further information and requests for resources, data, and codes should be directed to and will be fulfilled by the lead contact, Fumin Lei (leifm@ioz.ac.cn).

Materials availability

This study did not generate new unique reagents.

Experimental model and subject details section

Great tit and mice embryos were washed 30 times for 5 min with physiological saline. Internal organs were removed and sectioned with scissors. The pieces were then digested for 30 min with trypsin containing 0.25% EDTA (Thermo Fisher, 25200072). The supernatant was then centrifuged at 1,500 rpm for 5 min to harvest GEF and MEF. Cells were cultured in high glucose Dulbecco’s Modified Eagle’s medium (DMEM, Thermo Fisher, 10569010) supplemented with 10% fetal bovine serum (Thermo Fisher, 10099141) and a 1% penicillin/streptomycin mix (Thermo Fisher, 15140122) in an incubator at 37°C with a 5% CO2 atmosphere. The conditions of UVR exposure, low temperatures, and hypoxia treatments were chosen to simulate high-altitude environments based on previously published records. GEF and MEF were treated with UVR irradiation assays (Dahlback et al., 2007; Yu et al., 2016) comprising UVR irradiation at 100 J/m2 for 2 min using a UVR bench lamp with a medium wavelength of 308 nm. The cells were then harvested at 0 and 1 h after irradiation. GEF and MEF were subjected to low temperature assays (Beall, 2007; Ou et al., 2018) by placing them in a refrigerator at 4°C for 1 h or 4 h. Hypoxia pre-conditioning was performed by culturing cells in a tri-gas incubator for 4 h or 24 h with an oxygen concentration of 3% to simulate O2 levels at 2,900 m (Powell and Garcia, 2000). Untreated cell groups were used as controls.

Method details

Field collection of embryos was performed with the permission of the State Forestry Administration of China and conformed to the National Wildlife Conservation Law of China. All animals were maintained in accordance with guidelines established by the Institute of Zoology for the care and use of experimental animals. The experiments were also approved by the animal experimental and medical ethics committee of the Institute of Zoology at the Chinese Academy of Sciences.

RNA sequencing

The total RNA of cells was extracted using the TRIzol reagent (Invitrogen, USA). cDNA libraries were then prepared using standard Illumina protocols (Quail et al., 2009). Paired-end 150 bp libraries were sequenced with the Illumina HiSeq 4000 sequencing platforms. Primary quality control and secondary quality control were performed after generation of the sequencing FASTQ files. Primary quality control was conducted using our internal software. Default quality control settings were as including 1) removed of paired reads when the N content of single-ended sequencing reads exceeded 3% of the length of the read length, 2) removed of paired-end reads when the number of low-quality bases contained in a single-ended sequencing read exceeded 50% of the length of the read length, and 3) removed of reads when the adapter sequence did not match at least 8 bp of the expected sequence. Secondary quality control was conducted using a reference genome and the bwa program (Version: 0.7.15-r1140) (Jo and Koh, 2015). The first one million reads were extracted from the FASTQ file after primary quality control to compare against the input file. Mapp against reference rRNA sequences were conducted using Bowtie 2 (version 2.3.2) (Langmead and Salzberg, 2012) and the first ten thousand reads as the comparison input file. Gene expression levels were measured via reads per transcripts per kilobase of exon model per million mapped reads (TPM).

Weighted gene co-expression network analysis (WGCNA)

A WGCNA was performed using the WGCNA in R package (Langfelder and Horvath, 2008) to classify genes into modules. The analysis detects modules, or groups of correlated genes, which can reflect either shared regulation by common transcriptional factors or coordinated regulation by independent transcriptional factors with similar patterns of activation. Specifically, these groups were determined based on the short-term stress simulation experiments. The scale-free topologies of the networks were assessed for various β shrinkage parameter values according to the WGCNA user manual, wherein a β = 14 value provided a satisfactory fit to scale-free topology.

Transcriptional factor binding motifs were analyzed for each module using the CLOVER program (https://zlab.umassmed.edu/MotifViz/clover.html) to identify transcriptional factor binding site motifs overrepresented within the annotated proximal promoters of the genes within each module (Quach et al., 2016). The coding sequences of all genes in each module were used to analyze the common transcriptional factors for each module. The most abundant transcriptional factors for each module were subjected to a sequence LOGO analysis (http://jaspar.genereg.net) to evaluate sequence composition.

To evaluate putative gene functions, GO enrichment analysis and KEGG enrichment analyses were performed using g: Profiler (g: SCS threshold is the default parameter and p < 0.05 is the threshold. g: GOSt is the functional profiling.) (Raudvere et al., 2019) and Java GSEA4.0.3 (Reimand et al., 2019). Statistical significance of function enrichment was performed using Fisher’s Exact Test.

Differential expression analysis

Differential expression analysis was conducted by statistically comparing normalized read count data to determine quantitative changes in expression levels between treatment groups. DEGs were detected using a combination of DESeq2 (Love et al., 2014) packages. Genes with equal or above two-fold change in relative expression levels between the treatment and control groups, along with an FDR-adjusted p value < 0.05 were considered differentially expressed (Storey and Tibshirani, 2003). A MA plot was used to evaluate the differential expression results by plotting the log2 fold changes (y axis) against the mean expression signal (x axis) (Love et al., 2014). Volcano plots (Cui and Churchill, 2003) were also used to examine statistical significance (p values) versus fold change values.

Whole-genome resequencing and data processing

Genomic resequencing data of the great tit from highland and lowland habitats were downloaded from NCBI (PRJNA 274877) (Qu et al., 2015). Samples were obtained from a previous study of 32 great tits and were added 6 highland samples, including 12 individuals from the eastern Himalayas (highland), 16 from central/eastern China and 10 from inner Mongolia and the Mongolian lowlands. The data from the 38 great tits were mapped to the great tit genome using BWA (Li and Durbin, 2009). The alignment was processed using GATK to call SNPs (McKenna et al., 2010). SNPs were then filtered using VCFtools (Danecek et al., 2011) with parameters similar to our previous study (Qu et al., 2015). SNPs with a minor allele frequency (MAF) of less than 5% in the population were removed from the dataset.

Selection analyses

XP-CLR was used to detect signals of selective sweep, which involves jointly modeling the multilocus allele frequency between highland group and lowland group tits (Chen et al., 2010). XP-CLR scores were calculated using scripts available at http://genepath.med.harvard.edu/∼reich/ using the following parameters: sliding window size of 0.5 cM, a grid size of 2,000, maximum number of SNPs within a window size of 200 bp, and a correlation value for two SNPs weighted with a cutoff of 0.95. The highland group was considered as the object population, while the lowland group was considered the reference population. Haplotypes for each chromosome were inferred by plink (Purcell et al., 2007) and the general genetic map used for great tits was 1,917 cM = 1,200 Mbp (van Oers et al., 2014). With these, we first selected the regions with XP-CLR values in the top 5% of the empirical distribution. Next, the regions with the expected allele frequencies (Highland > Central East and Highland > Mongolian, or Highland

The XP-EHH (Sabeti et al., 2007) statistic was estimated for the highland and lowland groups. XP-EHH scores were calculated using scripts available at https://cran.r-project.org/web/packages/rehh. The highland group was again taken as the object population, while the lowland group was used as the reference population.

Positively selected genes were identified using the branch-site likelihood ratio test of nonsynonymous substitution rate to synonymous substitution rate (dN/dS) for each gene over 4 high-altitude vertebrates and their low-altitude relatives (Tang et al., 2017), 3 high-altitude passerine birds and their low-altitude relatives in PAML (Yang, 2007) from previous studies (Hao et al., 2019).

Evolutionary conservation analysis of MAPK1

The genes and pathway related to glucose metabolism of high-altitude adaptation in the great tit were retrieved from a previous study (Qu et al., 2015). The following pathways were retrieved: MAPK signaling pathway: Map 04010; PI3K-akt signaling pathway: Map 04151; mTOR signaling pathway: Map 04150; HIF-1 signaling pathway: Map 04066. PhastCons scores of the MAPK1 across 48 birds were extracted from Zhang (Zhang et al., 2014). MAPK1 sequences of the great tit were downloaded from NCBI. Sequence homology and phylogenetic analysis of amino acid sequences was conducted using MEGA. The phylogenetic tree was constructed based on MAPK1 from 14 different species. Domains in MAPK1 were identified via the Motif Scan (http://myhits.isb-sib.ch/cgi-bin/motif_scan) and NCBI conserved domains.

Protein 3D structure simulation

The 3D structures of the MAPK1 proteins of the great tit and humans were predicted using Phyre2 (Kelley et al., 2015), and then visualized using UCSF Chimera (Pettersen et al., 2004).

Identification of bird specific mutations

Protein sequences of MAPK1 were compared between human and birds. Two specific amino acid substitutions were defined as the human amino acid differed from that of all the other birds. We downloaded the MAPK1 protein sequence of humans and other birds from NCBI and compared them using Mega (Kumar et al., 2016). The conservative domains were identified by searching for conserved domains within coding nucleotide sequence (https://www.ncbi.nlm.nih.gov/Structure/cdd/wrpsb.cgi).

siRNA transfections

GEF were cultured in Dulbecco’s modified Eagle’s medium/Ham’s Nutrient Mixture F-12 supplemented with 10% fetal bovine serum (FBS) and 1% penicillin/streptomycin (Life Technologies, Carlsbad, CA, USA). The cells were plated at a density of 2 × 105 per mL in 24-well plates. For viability assays, 2 × 105 cells were seeded in 96-well plates and 25 mL cell culture flasks, respectively. The cells were maintained at 37°C in a 5% CO2 incubator and treated for the indicated intervals.

Transfections of MAPK1 siRNAs (Ruibo Co., Ltd, China) were performed using Lipofectamine 2000 (Invitrogen, 11668–027) according to manufacturer instructions. First, cells were plated the previous day and they were 70–90% confluent at the time of transfection, prepared plasmid siRNA-lipid complexes, and diluted suitable amounts of Lipofectamine Reagent and siRNA in Opti-MEM medium. Then, transfections were incubated for 5 min at room temperature and we added siRNA-lipid complexes to the cells. Twenty-four hours after transfection, cells were harvested for real-time quantitative reverse transcription PCR (RT-qPCR). Forty-eight hours after transfection, cells were harvested for western blotting analysis.

Cell apoptosis and cell cycle detection

The effects of siMAPK1 on cell apoptosis and cell cycle were evaluated using the Annexin V-FITC/PI Apoptosis Assay Kit (Vazyme Biotech Co., Ltd, NJ, USA) and the Cell Cycle Detection Kit (KeyGen Biotech, Nanjing, China), respectively, according to manufacturer protocols. Cells were resuspended in 1 × binding buffer to a final density of 1 × 105 to 5 × 105 cells/mL. Five microliters of Annexin V-FITC were added to each 100 μL of suspended cells, mixed and incubated for 20 min at room temperature in darkness. Then, 5 μL of propidium iodide staining solution was added, mixed and incubated at RT in darkness for 5 min. Finally, 400 μL of 1 × binding buffer was added, mixed and analyzed on a flow cytometer (FACSCalibur, BD, Franklin Lakes, NJ, USA). The same method was used to analyze cell cycle profiles. Three replicated experiments were performed.

Quantification and statistical analysis

All experimental data are presented as the means ± SEM of three independent experiments, and p < 0.05 was considered statistically significant. ANOVA and a subsequent Tukey’s post hoc test were performed to determine the significance level of each effect. All statistical analyses were performed using SPSS 17.0.

REAGENT or RESOURCE	SOURCE	IDENTIFIER
Chemicals, peptides, and recombinant proteins
trypsin	Thermo Fisher	25200072
DMEM	Thermo Fisher	10569010
fetal bovine serum	Thermo Fisher	10099141
penicillin/streptomycin	Thermo Fisher	15140122
Lipofectamine 2000	Invitrogen	11668–027
Annexin V-FITC/PI Apoptosis Assay Kit	Vazyme Biotech	A211-01
Cell Cycle Detection Kit	KeyGen	KGA
Deposited data
Raw data	This paper	https://www.ncbi.nlm.nih.gov/bioproject/PRJNA593713
Software and algorithms
GraphPad Prism	GraphPad	https://www.graphpad.com/scientific-software/prism/
Adobe Illustrator	Adobe	https://www.adobe.com/products/illustrator.html
Bowtie 2	Langmead and Salzberg (2012)	http://bowtie-bio.sourceforge.net/bowtie2/index.shtml
bwa	Jo and Koh (2015)	https://sourceforge.net/projects/bio-bwa/files/
WGCNA (R package)	Langfelder and Horvath (2008)	https://horvath.genetics.ucla.edu/html/CoexpressionNetwork/Rpackages/WGCNA/
g: Profiler	Raudvere et al. (2019)	https://biit.cs.ut.ee/gprofiler/gost
GSEA4.0.3	Reimand et al. (2019)	https://software.broadinstitute.org/cancer/software/gsea/wiki/index.php/GSEA_v4.0.x_Release_Notes
DESeq2	Love et al. (2014)	https://bioconductor.org/packages/release/bioc/html/DESeq2.html
GATK	McKenna et al. (2010)	https://gatk.broadinstitute.org/hc/en-us
VCFtools	Danecek et al. (2011)	http://vcftools.sourceforge.net/
XP-CLR	Chen et al. (2010)	https://github.com/hardingnj/xpclr
plink	Purcell et al. (2007)	https://www.softpedia.com/get/Science-CAD/PLINK.shtml
XP-EHH	Sabeti et al. (2007)	https://bio.tools/xp-ehh#!
PAML	Yang (2007)	http://abacus.gene.ucl.ac.uk/software/paml.html
Phyre2	Kelley et al. (2015)	http://www.sbg.bio.ic.ac.uk/phyre2/html/page.cgi?id=index
UCSF Chimera	Pettersen et al. (2004)	https://www.cgl.ucsf.edu/chimera/
Mega	Kumar et al. (2016)	https://www.megasoftware.net/dload_mac_beta