Comparative transcriptomic and metabolomic analysis reveals pectoralis highland adaptation across altitudinal songbirds

Abstract Pectoralis phenotypic variation plays a fundamental role in locomotion and thermogenesis in highland birds. However, its regulatory and metabolic mechanisms remain enigmatic to date. Here, we integrated phenomic, transcriptomic, and metabolomic approaches to determine muscle variation and its underpinning mechanisms across altitudinal songbirds. Phenomics confirmed that all highland birds had considerable increases in muscle oxidative capacity, capillarity, and mitochondrial abundance in our study. Correspondingly, transcriptomic analyses found that differentially expressed genes in phenotype‐associated modules enriched for blood vessel, muscle structure development, and mitochondrial organization. Despite similar traits and functional enrichments across highland birds, different mechanisms drove their occurrence in high‐altitude tree sparrow and 2 snow finches. Importantly, a metabolic feature shared by all the 3 highland birds is the improvement in insulin sensitivity and glucose utilization through activating insulin signaling pathway, which is vital to increase muscle oxidative capacity and maintain metabolic homeostasis. Nevertheless, fatty acid biosynthesis and oxidation are enhanced in only 2 snow finches which had a long evolutionary history on the high plateau, also differing from ketone body metabolism in recently introduced colonizer of the tree sparrow of the high plateau. Our study represents a vital contribution to reveal the regulatory and metabolic basis of pectoralis variation across altitudinal songbirds.


INTRODUCTION
At high altitude, hypoxia and hypothermia are 2 significantly severe challenges to aerobic exercise and thermogenesis in small endotherms, because a high rate of O 2 flux should be concurrently sustained to thermogenesis in cold temperature (Scott et al. 2015). Therefore, physiological modifications, including typically greater capillary density, more oxidative fiber, and higher proportion of subsarcolemmal mitochondria have been found in pectoralis of highland birds (Scott et al. 2009;Xiong et al. 2020). Understanding the molecular and biochemical mechanism of pectoralis modification is essential to investigate the functional evolution of phenotypic variation under hypoxic hypothermia. However, the regulatory and metabolic basis of muscle phenotype remains largely unknown for highland birds, in particular resident species with different evolutionary history.
In birds, pectoralis muscle is a fundamental tissue for locomotion and thermogenesis, while different fiber types (e.g. slow oxidative (SO), fast oxidative glycolytic (FOG), and fast glycolytic (FG) fiber) have been suggested to perform diverse functions due to their various contractile speed and metabolic capacity (Schiaffino & Reggiani 2011). Therefore, adaptive modifications in pectoralis of highland migrant birds are mainly reflected in changes in the relative proportion of fiber types (Scott et al. 2009). Previous studies, however, demonstrated that the pectoralis muscle of most small songbirds comprised exclusively FOG fibers (Marquez et al. 2006;Xiong et al. 2020). Thereby, small songbirds would take an alternative way to transform muscle phenotype but not to change muscle fiber types for responding to highland environmental stress. Moreover, resident small songbirds with different evolutionary time in the highland might be proposed to evolve various muscle characteristics.
Previous work mainly focused on uncovering bird genomic divergence to high altitude adaptation, while the role of the changes in gene expression contributing to muscle variation is given less attention (Qu et al. 2013(Qu et al. , 2020. In high-altitude rodents, the oxidative phenotype and capillarity of the gastrocnemius are related to the differential expression of genes involved in energy metabolism, muscle plasticity, and cell stress response (Scott et al. 2015). Hence, transcriptomic responses in pectoralis can prove an important insight into genes and biochemical pathways involved in trait variation or evolutionary adaptation to highland environment (Scott et al. 2015). Gene expression networks and key genes associated with phenotype are identified through weighted gene coexpression network analysis (WGCNA) in previous studies (Velotta et al. 2018;Hao et al. 2019). Likewise, skeletal muscle plays a central role on whole-body metabolism and also serves as an important contributor to maintain glucose homeostasis (Shepherd & Kahn 1999). Birds have higher plasma glucose concentrations than other vertebrates of similar body mass in part to the insulin insensitivity for the absence of Glut4 which can de-crease blood glucose after feeding in mammals (Polakof et al. 2011;Xiong & Lei 2021). Therefore, how to maintain energetic and redox homeostasis in pectoralis would be a great challenge for birds surviving under hypoxic hypothermia.
Snow finches (Onychostruthus taczanowskii and Pyrgilauda ruficollis) and tree sparrow (Passer montanus) are closely related within the Old World sparrows (Passeridae) (Päckert et al. 2021). Snow finches are endemic species on the Qinghai−Tibet Plateau (QTP) with a strict elevation distribution range from 3500 to 5100 m, whereas the tree sparrow is an introduced colonizer on the QTP with altitudinal distribution from sea level to 4400 m, which migrated with Tibetans as previous study (Fig. 1a) (del Hoyo et al. 2018;Qu et al. 2020). Snow finches and high-altitude tree sparrow separately have relative longer to shorter evolutionary time to the QTP, which comprise a powerful model system to determine the pectoralis variation and its underlying genetic and metabolic basis under the high-altitude environments.
Based on this model system, we combined phenomic, transcriptomic, and metabolomics measures to the following 3 objectives in this study. First, we found evidences of muscle differences on fiber composition, capillarity, mitochondrial distribution, and abundance across altitudinal songbirds through histochemical methods. Second, we utilized a functional transcriptome approach to decide the regulatory basis of these variations through the differentially expressed genes (DEGs) and co-expression network analysis (WGCNA). Third, we used highly sensitive gas chromatography-mass spectrometry (GC-MS) combined with the RNA-seq data to identify their own metabolic mechanisms.

Sampling
Wild birds were caught with mist-net in summer of 2016, highland bird populations from Qinghai Province about 3200 m (37.033°N, 99.73°E, 10 tree sparrows, P. montanus) and 3900 m (35.4°N, 99.25°E, 10 rufousnecked snow finch, P. ruficollis and 10 white-rumped snow finch, O. taczanowskii) and lowland population (40.39°N, 116.67°E, n = 8 for transcriptome and 10 for GC-MS analysis) of the tree sparrow from Yanqi Lake of Beijing at 80 m (Table S1, Supporting Information). 100-200 μL of bird blood was collected from wing vein, centrifuged 10 min at 2500 rpm and separately stored at −80°C. Birds were euthanized Sampling elevations and body mass of 3 bird species, histological analysis, and electron microscope structure of pectoralis major. (a) Sampling elevations and body mass of 3 species. Phylogenetic tree is reconstructed based on complete mitochondrial genome. Maps were generated using ArcGis (https://www.esri.com). (b) Oxidative fibers and capillaries were identified via staining in succinate dehydrogenase activity and alkaline phosphatase activity, respectively; representative transmission electron micrographs for pectoralis major. White arrow, muscle fiber; black arrow, capillary; white arrowhead, mitochondria; black arrowhead, myofibril. through heart-compressed method after blood collection when sampling. Pectoral major muscle was dissected immediately following euthanasia, flash-frozen in liquid nitrogen, and stored at −80°C. All procedures performed on animal were approved by Institute of Zoology Animal Care Committee.

Muscle histology and transmission electron microscopy
Oxidative muscle type and capillarity were performed following Scott and Johnston (2012). Briefly, pectoral major muscle was dissected and samples were taken third way along the length of the sternum, covered with OCT, and frozen in isopentane (cooled in liquid N 2 ). Muscle sections (10 mm) were obtained in a Cryostat Microtome (Leica CM900, Germany) at −20 • C. More details were presented in the Supporting Information.
The flight muscle was removed from an intermediate depth and then fixed at 4°C for 24-48 h in 2% glutaraldehyde in 0.1 M PBS buffer at pH 7.4. Small muscle blocks (2 mm × 2 mm) were prepared and post-fixed in 1% osmium tetroxide in 0.1 M PBS buffer for 1 h, dehydrated in a graded ethanol series (50%, 70%, 70%, 95%, 95%, 100%, 100%), and embedded in epoxy resin. Ultra-thin sections were cut on a Leica UC7 ultramicrotome and placed on copper grids. The sections were post-stained with uranyl acetate and lead citrate. Images were collected using a transmission electron microscope (Tecnai G2 F20 TWIN TMP, USA). We measured mitochondrial volume density using stereological methods as previously described (Weibel 1979). Grid size of 90 nm was used to estimate mitochondria and lipid droplet volume at a square size of 4460 × 4460 nm.

Plasma glucose, lactate, insulin, tissue glycogen measurement, and enzyme activity assays
To quantify circulating glucose levels, fasting plasma was obtained from 8 birds from each species and was measured using an Accu-Chek blood glucose meter (Roche Diagnostics). Insulin content of fasting plasma was performed using insulin ELISA kit (Cloud-Clone), lactate in plasma and muscle glycogen were determined with a relative assay kit (Solarbio) according to the manufacturer's instructions. More information about enzyme activity is presented in Supporting Information.

RNA isolation, library preparation, sequencing, quantification, and normalization
Total RNA was extracted from each pectroralis major muscle of the 38 samples (8 lowland tree sparrows, 10 highland tree sparrows, 10 rufous-necked snow finches, and 10 white-rumped snow finches) using Trizol RNA isolation reagents (Invitrogen Corp., Carlsbad, CA). RNA integrity was assessed using the RNA Nano 6000 Assay Kit of the Agilent Bioanalyzer 2100 system (Agilent Technologies, CA, USA). A total amount of 3 μg RNA per sample was used as input material for the RNA sample preparations. Sequencing libraries were generated using NEBNext UltraTM RNA Library Prep Kit for Illumina (NEB, USA) following manufacturer's recommendations and index codes were added to attribute sequences to each sample.
The clustering of the index-coded samples was performed on a cBot Cluster Generation System using TruSeq PE Cluster Kit v3-cBot-HS (Illumia) according to the manufacturer's instructions. After cluster generation, the library preparations were sequenced on an Illumina Hiseq platform and 125 bp/150 bp paired-end reads were generated. Trimmomatic (Bolger et al. 2014) was used to filter reads containing adapter, reads containing ploy-N and low quality reads based on read quality checked with FASTQC (Andrews 2016). The parameters used were as follows: sliding window = 4-bp; Phred 33 quality scores = 20; Min read length = 50. Adapter sequences, if detected, were removed. Clean data with high quality was mapped from each species to respective genomes using STAR with default parameters (Dobin et al. 2013). We used the reciprocal best-hit method to generate tree sparrow-rufous-necked snowfinch and tree sparrow-white-rumped snowfinch orthologs, respectively. The orthologs shared by 3 species were obtained by intersecting the lists of above 2 orthologs. After the reads were mapped to the reference genomes (Qu et al. 2020), expression quantifications of genes and transcripts were performed using RSEM (Li & Dewey 2011). Expression levels for genes with one-to-one orthologs in all 3 bird species (n = 12 951) were normalized with a RLE (relative log expression) method across muscle samples (Maza 2016).

Differential expression analysis and weighted gene coexpression network analysis
To identify genes related with pectoralis variation across altitudinal songbirds, we performed differential expression analysis. Differentially expressed genes (DEGs) were calculated based on the Negative Binomial distribution and independent filtering was enabled in a R/Bioconductor package Deseq 2 (R version 3.5.1) with a false discovery rate (FDR) < 0.05 based on Benjamini-Hochberg method to control the FDR in multiple tests context in identifying significantly differentially expressed genes (Love et al. 2014). The cut off values for log-2-fold change were set at 0.59 and −0.59. Only genes with count number greater than 1 in at least 4 samples were included in differential gene expression analysis.
Then we performed a weighted gene co-expression network analysis via WGCNA to identify gene modules associated with muscle phenotypic variation and its potentially genetic basis. Only inter-and intra-species differentially expressed genes (DEGs) were included in the analysis. Briefly, we constructed the weighted gene co-expression network using the normalized, log 2transformed counts to analyze the DEGs with the block-wiseModules function in WGCNA (Langfelder et al. 2008) for inter-and intra-species, respectively. An adjacency correlation matrix is calculated for the DEGs, and the correlations are weighted to a soft threshold power β which favors strong correlations over weak one (Zhang & Horvath 2005). For each pair of genes, a robust measure of network interconnectedness is calculated based on the adjacency matrix. For our analysis, the parameters used were as follows: for population network, maximum block size = 966 genes, power (β) = 18; for species network, maximum block size = 2457 genes, power (β) = 22; minimum module size = 25; minimum height for merging modules = 0.25; maximum height for cutting the tree = 0.90. The remaining parameters were kept at the default settings.
Co-expression modules associated with phenotypic gradients were identified using a principal component analysis (PCA) of gene expression with the blockwise-Modules function in WGCNA. Each module was summarized by an eigengene, which is the first principal component of the scaled module expression. Thus, the module eigengene explained the maximum amount of variation of the module expression levels. A Student's asymptotic test with the corPvalueStudent function was used to determine P values of the correlation. To specify genes potentially explaining muscle phenotypic variation, we identified hub genes which may be central to the architecture of the regulatory networks represented by each co-expression module. The hub genes are calculated by their first principal component (PC1), the module eigengene (a summary of overall module expression) (Han et al. 2004). We then implemented gene ontology categories within modules pos-itively and negatively correlated with muscle traits using G:PROFILER (FDR < 0.05) (Reimand et al. 2007). According to the results of functional GO enrichment analysis, intramodular hub genes were identified as candidate genes driving potentially muscle variation.

Statistical analyses
The 2-tailed student's t-test was used to determine statistical significance. Correlative analysis was carried out using linear regression. For all tables and figures, P value < 0.05 was considered significant.

Pectoralis phenotypes and PCA
Pectoralis major and body mass of all highland birds were significantly larger than those of lowland tree sparrow, while only highland tree sparrow had a bigger relative ratio of pectoralis major mass to body mass ( Fig. 1a and Table 1), reflecting pectoralis as a contributor to the increase of body mass of them. Histological analysis found that highlanders evolved a higher oxidative capacity of pectoralis through increasing the sizes of fiber for highland tree sparrow or myofibril for 2 snow finches ( Fig. 1b and Table 1). In contrast to the increased musclefiber size in highland tree sparrows, snow finches had thicker myofibril (Table 1), indicating that different modifications of muscle phenotype occur in intraspecies (highaltitude tree sparrow vs. low-altitude tree sparrow) and interspecies (2 snow finches vs. low-altitude tree sparrow). Additionally, the ratio of wing length to body length significantly associated with muscle fiber area (R 2 = 0.36, P = 0.0002, Fig. S6a, Supporting Information) and myofibril diameter (R 2 = 0.51, P < 0.0001, Fig. S6b, Supporting Information).
The improvement in oxidative capacity of flight muscle is also attributed to more capillarity and mitochondrial volume density in highland birds. All highland birds exhibited an enhanced capacity of oxygen transport though the higher capillary to fiber ratio (CF) and capillary size ( Fig. 1b and Table 1). Snow finches also had a significantly great capillary density (CD) in comparing to lowland tree sparrows (Table 1). Highland tree sparrow did not increase CD due to a decreased fiber density (Table 1), indicating that increased CF is a compensatory modification for diminishing of fiber density.
Additionally, highland birds had a greater proportion of subsarcolemmal mitochondria (DS) and a greater total mitochondrial volume density (Vv (mt)) in pectoralis ( Fig. 1b and Table 1). The last trait in the tree sparrow was entirely contributed to the greater volume density of subsarcolemmal mitochondria (Vv (ssm)), while snow finches had the greater Vv (ssm) and volume density of intermyofibrillar mitochondria (Vv (imm)) ( Table 1). Moreover, the Vv (mt) was connected positively with the Vv (LD) (R 2 = 0.44, P = 0.0004, Fig. S6c, Supporting Information) and myofibril size (R 2 = 0.56, P < 0.0001, Fig. S6d, Supporting Information). The difference in mitochondrial abundance between highland and lowland species was greater than the difference between altitudinal populations of the tree sparrow, reflecting that oxidative capacity of flight muscle was greater in snow finches than highland tree sparrows.
The first 3 principal components of the PCA on the basis of 12 phenotypic traits explained most of variations in pectoralis across inter-and intra-species (Fig. S1, Supporting Information). These traits can be grouped into 5 classes including body weight, lipid volume density, muscle fiber or myofibril (FD, FS and MD), capillarity (CD, CA and CF), and mitochondria (Vv (mt), Vv (ssm), Vv (imm), and Vv (sim)). Within species, the PC1 (59.1%) represented mainly most of the variance in mitochondrial volume density, lipid droplet, capillarity, and weight. PC2 (13.1%) and PC3 (8%) were also explained by myofibril and volume density of intermyofibrillar mitochondria. In addition, PC1 (58.7%) among species showed variance of capillarity, myofibril, and mitochondrial morphology, while PC2 (18.4%) and PC3 (6.5%) explained muscle fiber density, capillary size, as well as lipid droplet, respectively. These results showed that the PC1 of PCA could reflect most of muscle phenotype variations.

Gene expression differences and WGCNA
We sequenced 38 samples of pectroralis major muscle for 2 snow finches and 2 tree sparrow populations, which had the matching muscle phenotype. The G20 and G30 of each species and populations were more than 98% and 90%, and GC content was all around 49-50%. After quality control, we generated a total of 141 Gb of high-quality RNA-seq data for subsequent analyses. Samples were aligned to their respective reference genomes, with overall mapping rates ranging from 87.51% to 94.58% across 38 samples. The mean alignment rates were 90.83%, 90.98%, 92.55%, and 91.92% for lowland tree sparrow, highland tree sparrow, rufous-necked snow finch, and white-rumped snow finch, respectively. Therefore, our sequencing data was of high quality and was robust for differentially expressed analyses (Table S1, Supporting Information).
Weighted gene correlation network analysis using WGCNA organized 966 DEGs within populations and 2457 shared DEGs in snow finches to 4 and 5 modules, respectively (Fig. 2a,b; Figs S2-S4, Supporting Information). In population network (high-and low-altitude tree sparrows), blue module and brown module correlated positively with PC1 (Fig. 2c). The DEGs in blue module were enriched into "muscle structure development" (P < 0.001) and "muscle contraction" (P = 0.0281) (Table S4, Supporting Information). The DEGs in brown module were enriched into "mitochondrion" (P < 0.001), "respiratory chain" (P < 0.001) and "oxidative phosphorylation" pathway (P = 0.001) (Table S4, Supporting Information). Also, yellow module and turquoise module associated negatively with PC1 (Fig. 2b). The DEGs in these modules were separately enriched into "organelle subcompartment" (P = 0.0300) and "muscle system process" (P = 0.001), "mitochondrion" (P = 0.023) as well as some metabolic pathways (e.g. "carbohydrate derivative metabolic process" [P < 0.001] and "citrate cycle" [P = 0.037]) (Table S4, Supporting Information). As to species network, blue module and turquoise module associated positively with PC1 (Fig. 2d). The DEGs in turquoise module were enriched into many developmental processes including "muscle structure development" (P < 0.001) and "vasculature development" (P < 0.001), "mitochondrion" (P < 0.001), and some metabolic processes such as "lipid metabolic process" (P = 0.009) and Figure 2 Co-expression analyses of differentially expressed genes and muscle traits via WGCNA. Average linkage clustering tree and correlation between co-expression modules and muscle traits across intra-species (a) as well as inter-species (b). (c) Co-expression modules identified via WGCNA. Four modules are associated significantly with PC1 of pectoralis major variation between populations (2 positively connected modules, blue module: 300 genes, brown module: 50 genes; 2 negatively connected modules, yellow module: 51 genes, turquoise module: 438 genes). (d) Four modules are correlated with PC1 of pectoralis major variation across species (2 positively associated module, turquoise module: 1417 genes, blue module: 181 genes; 2 negatively connected modules, brown module: 113 genes, yellow module: 103 genes). Green module (68 genes) is associated negatively with PC2. Network view of angiogenesis, muscle development, and mitochondrial organization in correlated modules in intra-species (e) and inter-species (f) level. Shapes indicate modules and are colored according to GO terms. Shape size is presented in accordance with module membership, and the line border color reflects the correlation with PC1. The co-expression between genes is reflected by edge line width. "inositol phosphate metabolism" (P = 0.001) (Table S4, Supporting Information), and yellow module and brown module associated negatively with PC1 (Fig. 2d). The DEGs in these modules were respectively enriched into "organelle" (P = 0.009) and "membrane-bounded organelle" (P < 0.001) (Table S4, Supporting Information).
Only green module had a negative correlation with PC2 (Fig. 2b) and genes in this module were related to "membrane part" (P = 0.004) ( Table S4, Supporting Information). These gene ontology (GO) enrichment and KEGG pathways associated with angiogenesis, muscle development, mitochondrial organization, and metabolic process, suggested the potential regulatory basis of muscle phenotypes.
As to species network, 2 hub genes of the blue module were involved in vasculature development (MEF2C), and muscle structure development (MEF2C, ASB2) (Figs 2f,3c; Table S5, Supporting Information). The turquoise module had 12 hub genes, 11 of these 1170 (expect for WIPI2) only in snow finches but not in highland tree sparrow comparing with lowlander, probably indicating the regulatory difference of pectoralis phenotypes between snow finches and highland tree sparrow. The 5 genes in this 12 hub gene subset were involved in angiogenesis (EPAS1, VEGFR3, VEGFR4, PIK3C2A, KIAA1462), 3 were involved in muscle structure development (FLNB, MAML1, EPAS1), and 6 were involved in mitochondrion organization (LIG3, HIGD2A, WIPI2, EPAS1, MID51, HTT) (Figs 2f,3c; Table S5, Supporting Information). EPAS1, comprising a dimeric complex with ARNT responding to hypoxia, was upregulated in snow finches and played a central role on development process including angiogenesis, muscle development, and mitochondrion organization in species network ( Fig. 2f; Table S6, Supporting Information) (Rankin et al. 2008;Skuli et al. 2009;Park et al. 2016). In addition, several key downstream genes of EPSA1 including VEGFR3, and VEGFR4, mediating angiogenesis as well blood vessel morphogenesis, were also overexpressed hub genes and associated with CD (VEGFR3, R 2 = 0.28, P = 0.0074) and CF (VEGFR4, R 2 = 0.62, P < 0.0001) across species ( Fig. 3c; Table S6, Supporting Information). MID51 involving in mitochondrial fission was a hub gene with a higher degree and its expression positively associated with the Vv (mt) (R 2 = 0.74, P < 0.0001). In these modules, some genes would have not high intramodular connectivity and degree via WGCNA, while function enrichment analyses and association test with muscle phenotypes also revealed an essential impact on regulating skeletal muscle development, angiogenesis, and mitochondrial biogenesis (Figs 2f,3a; Table S6, Supporting Information).
To confirm the real expression of some hub DEGs from protein level, we performed western blots for 2 key genes including MEFC and EPAS1, for an example. In accord with normalized transcript counts, protein levels of MEF2C were improved in all highland songbirds compared to lowlanders, while EPAS1 only in rufous-necked snow finch and white-rumped snow finch but not highland tree sparrows (Fig. S7, Supporting Information) ).

Metabolic basis of phenotypic variations
Given that skeletal muscle plays a central role on whole-body metabolism, we used highly sensitive GC-MS combined with the RNA-seq data to identify metabolic features in different muscle phenotypes. Fasting plasma glucose was lower in most of highland birds, while lots of monosaccharides and glycogen were accumulated in muscle fiber (Fig. 4a,b; Fig. S5a, Supporting Information), indicating improvement in glucose uptake and utilization. Highland songbirds including whiterumped snow finch and tree sparrow have significant low plasma lactate despite rufous-necked snow finch slightly decrease lactate (Fig. 4a). Although plasma insulin concentration was similar to those of lowland birds (Fig. 4a), highland birds had some DEGs involved in regulating insulin sensitivity and insulin signaling pathway (AR/ERα, PIK3C2A, PIK3C2B, SLC2A12 [Glut12], and PTP1B) (Fig. 4c). Additionally, expression and activities of some glycolytic enzymes (hexokinase, pyruvate kinase, and lactate dehydrogenase) were increased in highland birds (Fig. 4a,c). Consistently, glycolytic intermediates increased in all highland birds, with increased glucose and fructose-6-phosphate (Fig. 4b). Activity of citrate synthase (CS) was higher in highland birds (Fig. 4a), which supported a great Vv (mt) in pectoralis (Table 1). Meanwhile, concentrations of 6-, 5-, and 4-carbon intermediates downstream of CS (aconitate, α-ketoglutarate, fumarate, and malate) were increased (Fig. 4b). Collectively, the improvement of insulin sensitivity and glucose utilization was a shared feature of muscle metabolism to hypoxia through activating insulin signaling pathway.
Muscle phenotype suggested that snow finches had greater capacities of oxygen transport and mitochondrial oxidation, which probably predicted alternative fuel substrates in snowfinches. Consistently, muscle metabolite levels also indicated that total fatty acid concentration increased in 2 snow finches ( Fig. 5a; Fig. S5b, Supporting Information). The ratio of short-chain fatty acid (SFA), however, was great in snow finches instead of a low ratio of desaturated fatty acid (DFA) (Fig. 5a). Meanwhile, the ratio of long-chain fatty acid (LFA) also significantly decreased in snow finches (Fig. 5a), indicating DFA and LFA might be important fuel substrates for thermogenesis under low temperature. Additionally, PPARδ had a low expression level in snow finches, while DGAT2 activating lipid depot formation was up-regulated and its expression was positively associated with volume density of lipid droplet (R 2 = 0.49, P = 0.0002) ( Fig. 5b; Table  S6, Supporting Information). Interestingly, PGC1a and some PPAR target genes, involving in fatty acid transport, ω-oxidation, and lipid droplet degradation, had an increased expression and activity of HADH was also dominantly increased (Fig. 5a-c). Together, snow finches likely decrease PPARδ expression to promote intramuscular lipid biosynthesis and increase mitochondrial metabolism for fatty acid oxidation for suffering a long highland life history. In contrast, hydroxybutyrate (HB) decreased significantly in highland tree sparrow (Fig. 5a), indicating that ketone instead of fatty acid was a major anaplerotic substrate for energetic supply in a short evolutionary time to hypoxia.

DISCUSSION
By integrating a comprehensive array of analytical techniques, we illuminate muscle phenotypic variation and its underpinning regulatory as well as metabolic 1172 Figure 5 Fatty acid biosynthesis and oxidation in pectoralis. (a) 3-hydroxyacyl-CoA dehydrogenase (Hadh) activity, hydroxybutyric acid content, total fatty acid content, long chain/total fatty acid ratio, short chain/total fatty acid ratio, and unsaturated/total fatty acid ratio in pectoralis. * significant increase; † significant decrease. (b) Expression of genes for fatty acid transport, ω-oxidation, fatty acid degradation, and regulation of fatty acid metabolism were differentially expressed in snow finches. (c) Schematic summary of fatty acid metabolism in pectoralis. Black arrow, lipid droplet; black arrowhead, mitochondria. mechanism across altitudinal songbirds with different highland lifetime history. Despite considerable increases in muscle oxidative capacity, capillarity and mitochondrial abundance when facing hypothermia and hypoxia, different regulatory mechanisms might drive their occurrence between snow finches and high-altitude tree sparrow. A metabolic feature shared by all highland birds is the improvement in insulin sensitivity and glucose utilization through activating insulin signaling pathway, which is vital to maintain metabolic homeostasis in highland birds. Additionally, an alternative fuel would be taken by highland birds with different lifetime under hypothermia environment.

Pectoralis variations for highland survival
Great body mass and muscle weight in highland animals are in concordance with Bergmann's rule (Sun et al. 2017;Fan et al. 2019). Avian pectoralis is used to support the high metabolic costs through aerobic metabolism for flight or thermogenesis (Marquez et al. 2006). Highaltitude birds, however, suffer from severely limited oxy-gen supplement due to decrease in the partial pressure of oxygen. Therefore, the physiological responses to hypoxic hypothermia of highland are generally found by increasing oxidative capacity as well as oxygen delivery and by altering oxygen utilization in pectoralis (Scott et al. 2009;Xiong et al. 2020). Either size modification of muscle fiber or myofibril in highland birds plays an important role in heat production through shivering thermogenesis for maintenance of a constant core body temperature (Chaillou 2018). Additionally, various modifications of muscle phenotype across altitudinal birds are probably associated with the flight performance in hypoxic environment. For example, an increased size of myofibril would produce more contractile force to prevent flightless, probably accounting for why highland tree sparrows perform a weak flight capacity (Sun et al. 2016).
The improvement in oxidative capacity of pectoralis also attributes to more capillarity and mitochondrial content in highland birds. The increases in CD and CF are thought to enhance the diffusion capacity of oxygen from blood to muscle (Scott & Milsom 2007;Xiong et al. 2020). Capillary size increase may reduce capillary flow resistance directly and maintain capillary rheology in hypoxic environment (Hauck et al. 2004). Mitochondrion is the main place of aerobic oxidation and takes a majority of ATP production in eukaryote (Rich 2003). An increase in DS in highland birds may enhance ability to produce proton-motive force as reported in previous studies (Mahalingam et al. 2017). Overall, capacities of oxidation and oxygen delivery in pectoralis are prompted through differentially phenotypic variations in highland tree sparrow and snow finches, probably resorting from different adaptation time to hypoxic hypothermia.

Potential basis of pectoralis variations in differential gene expression
Given that the differences in muscle fiber and myofibril occur in highland tree sparrow and snow finches comparing to lowland tree sparrows, respectively, we propose the potential regulatory basis of pectoralis variation through DEGs and WGCNA analyses. Some genes are found to be involved in the regulation of muscle development via MRF4-MEF2 axis, which is thought to be a novel pathway regulating muscle fiber size and muscle mass in adult muscle tissue (Moretti et al. 2016). MRF4 and MEF2C are known to play a role in myogenesis (Hernandez-Hernandez et al. 2017). MRF4 knockdown in skeletal muscle of adult rodent stimulates MEF2 transcriptional activity and then causes muscle hypertrophy (Moretti et al. 2016). However, HDAC3, which can deacetylate and repress MEF2C, can be up-regulated by hypoxia-induced factor (HIF-1α), and increases cardiomyocyte hyperplasia in mice (Wu et al. 2011). Therefore, MRF4-MEF2 axis likely controls muscle hypertrophy in highland tree sparrow and causes muscle hyperplasia under overexpression of HDAC3 induced by EPAS1 in snow finches, which may also be a regulative process of postnatal but not embryonic myogenesis (Moretti et al. 2016).
All highland birds have an increase in capillary size along with the overexpression of SEMA3C which is also a hub gene in the highland tree sparrow. Therefore, SEMA3C probably increases integrin activation to control blood vessel size in highland birds as previous study proposed (Banu et al. 2006). Additionally, snow finches also have an increase in capillary density along with the overexpression of EPAS1 and its downstream genes. Mounting evidence suggests that EPAS1 has an important role in angiogenesis, though most of the studies are focused on HIF1α (Rankin et al. 2008;Skuli et al. 2009;Park et al. 2016). Overexpression of EPAS1 and ARNT in pectoralis of snow finches probably maintains an enhanced capillarity via controlling expression of target genes both in VEGF/VEGFR and ANG/TIE pathways (Rankin et al. 2008;Skuli et al. 2009;Park et al. 2016). Importantly, the increased expression of VEGFR3 promotes the expression of DLL4 and HEY2, and then controls the conversion of tip cell to stalk cell by reinforcing Notch signaling during angiogenic sprouting (Hellstrom et al. 2007). Nevertheless, the up-regulated expression of TIE2 and ANGPTL1 stimulates basement membrane deposition and pericyte detachment, thereby mediating vessel maturation (Augustin et al. 2009).
Autophagy is a lysosome-mediated degradation process for cytoplasmic components, and also plays an important role in mitochondrial homeostasis (Lee et al. 2012). Highland birds decrease the expression of genes involved in mitophagy, including WIPI2, LC3, and BNIP3 (Fig. 3a), indicating a common central role in maintaining mitochondrial volume density. Also, TRAK2 and CLUH are shared regulation genes of mitochondrial distribution in highland birds (MacAskill & Kittler 2010). However, differentially expressed genes for mitochondrial fission and fusion are various. Snow finches and highland tree sparrow induce mitochondrial fission via MID51 and MID49, respectively (Ni et al. 2015). Additionally, MIGA1 and MIGA2 also enhance mitochondrial fusion only in snow finches (Zhang et al. 2016). Besides, KLF4/ERR/PGC1a transcriptional complex regulates a broad spectrum of genes involved in mitochondrial biogenesis, dynamics, and metabolism (Dorn et al. 2015). All 3 genes are over-expressed in snow finches, but only ERR has an improved expression in the highland tree sparrow, probably indicating an increase in Vv (ssm).

Insulin sensitivity, glucose utilization, and fatty acid metabolism
One of the prevailing views considers energetic needs for thermogenesis in highland coming from glucose oxidation which requires less O 2 than fatty acid oxidation (Schippers et al. 2012). Glucose uptake from the bloodstream into cells is mediated by a family of facilitative glucose transporters (Gluts) (Shepherd & Kahn 1999;Xiong & Lei 2021). It has been long suggested that a low sensitivity to insulin is due to the absence of 1174 Glut4 in birds (Polakof et al. 2011). However, recent studies elucidate that Glut12 overexpression parallels that of Glut4, improves insulin sensitivity, and enhances glucose uptake in skeletal muscles of rodent and chicken (Coudert et al. 2015;Xiong & Lei 2021). Highland birds exhibit a similar concentration of plasma insulin and an accelerated glucose uptake in flight muscle along with high expression of Glut12, indicating an improvement of insulin sensitivity. Low plasma lactate in highland songbirds was likely driven by high expressions of lactate dehydrogenase (LDH) and also by high activity of LDH in highland tree sparrow. Previous studies have suggested that AR/ERα, PIK3C2A, PIK3C2B, and PTP1B mediate insulin sensitivity through activating insulin signaling pathway (Yu et al. 2013). It is the first time to find an enhanced glucose metabolism via activating insulin signaling pathway by gene expression in highland vertebrates. Additionally, increase in insulin sensitivity improves glucose oxidation (Yan 2018), which mainly exhibits in increased capacities of glycolysis and TCA cycle (Horscroft et al. 2017).
Snow finches also improve the capacities of fatty acid biosynthesis and oxidation, differing from ketone metabolism of the highland tree sparrow. Per mole of fatty acid oxidation can yield more ATPs than glucose despite more oxygen consumption (Schippers et al. 2012); thereby, it is an important energy strategy taken by small animals to cold exposure (Cheviron et al. 2012). The enhanced capacities of lipid oxidation probably occurs in snow finches facing prolonged cold and originates from the improvement of β-oxidation as well as ω-oxidation through overexpression of target genes or higher enzyme activity, which is likely regulated by PGC1a (Cheviron et al. 2012;Horscroft et al. 2017).

CONCLUSIONS
Determining the molecular and functional mechanisms underlying phenotypic variation is a fundamental goal of evolutionary biology and ecology. In this study, we confirm that highland birds have an extensive increase in capillarity, mitochondrial abundance, and oxidative capacity of muscle fiber. RNA-seq analyses reveal differentially expressed genes in modules associated with phenotypes enriched in blood vessel, muscle structure development, and mitochondrial organization. Despite similar phenotypes and functional enrichments across highland birds, different molecular and biochemical mechanism drive their occurrence in part due to various evolutionary histories. Moreover, a shared metabolic feature among high-land birds is the improvement in insulin sensitivity and glucose utilization through activating insulin signaling pathway, which is vital to maintaining metabolic homeostasis in highland birds. However, fatty acid biosynthesis and oxidation are enhanced in snow finches with a long evolutionary history, also differing from ketone body metabolism in recently introduced colonizer. To our knowledge, this study is the first to investigate pectoralis variation and its regulatory as well as metabolic mechanisms across altitudinal songbirds. Zhang YP, Liu XM, Bai J et al. (2016). Mitoguardin regulates mitochondrial fusion through MitoPLD and is required for neuronal homeostasis. Molecular Cell 61, 111-24.

SUPPLEMENTARY MATERIALS
Additional supporting information may be found online in the Supporting Information section at the end of the article. Figure S1 Rotated principal component analysis of flight muscle. Figure S2 18 and 22 of soft power for scalefree network construction of intraspecies (a) and interspecies (b) respectively. Figure S3 Modules correspond to muscle samples within species (a) and between species (b). Figure S4 Network heatmaps of the 966 and 2457 genes differentially expressed between intraspecies (a) and interspecies (b) respectively.

Figure S5
Relative levels of intermediates and fatty acids in fasting plasma. * significant increase; ‡ significant decrease. Figure S6 Associations between wing length/body length and fiber area (a), between wing length/body length and myofibril diameter (b), between Vv (mt) and Vv (LD) (c), between Vv (mt) and myofibril diameter (d), and between Vv (LD) and myofibril diameter (e). Figure S7 Protein levels of MEF2C and EPAS1. Table S1 Specimen information used for phenomic, transcriptomic and metabolomic analyses Table S2  Table S3 Gene ontology terms, Human Phenotype Ontology, and KEEG pathways enriched in shared differentially expressed genes Table S4 Gene ontology terms, Human Phenotype Ontology, and KEEG pathways enriched in muscle functional module of differentially expressed genes