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

doi:10.21203/rs.2.22368/v2

Download PDF

Research article

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

https://doi.org/10.21203/rs.2.22368/v2

This work is licensed under a CC BY 4.0 License

Journal Publication

published 21 Dec, 2021

Read the published version in Integrative Zoology →

Version 2

posted

You are reading this latest preprint version

Background

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.

Results

Phenomics revealed that all highland birds had considerable increases in muscle oxidative capacity, capillarity, and mitochondrial abundance. Correspondingly, transcriptomic analyses found that differentially expressed genes in modules associated with phenotypes enriched in blood vessel, muscle structure development, and mitochondrial organization. Despite similar traits and functional enrichments across highland birds, different mechanisms drove their occurrence in part for their own various evolutionary histories. Importantly, a metabolic feature shared by 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 species with a long evolutionary history, also differing from ketone body metabolism in recently introduced colonizer.

Conclusions

Our study represents a vital contribution to reveal the regulatory and metabolic basis of pectoralis variation across altitudinal songbirds.

Epigenetics & Genomics

pectoralis

hypoxia

altitude

gas chromatography-mass spectrometry

RNA-seq

At high altitude, hypoxia and hypothermia are two significantly severe challenges to aerobic exercise and thermogenesis in small endotherms, because a high rate of O₂ flux should be concurrently sustained to thermogenesis in cold temperature[1]. 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 [2]. 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 [3]. Therefore, adaptive modifications in pectoralis of highland migrant birds are mainly reflected in changes in the relative proportion of fiber types [2]. Previous studies, however, demonstrated that the pectoralis muscle of most small songbirds comprised exclusively of FOG fibers [4]. Thereby small songbirds would take an alternative way to transform muscle phenotype 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, whilst the role of the changes in gene expression contributing to muscle variation is given less attention [5]. 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 [1]. 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 [1]. Likewise, skeletal muscle plays a central role on whole-body metabolism and also serves as an important contributor to maintain glucose homeostasis [6]. Birds have higher plasma glucose concentrations than other vertebrates of similar body mass in part to the insulin insensitivity for the absence of Glut 4 [7]. Therefore, how to maintain energetic and redox homeostasis in pectoralis would be a great challenge for birds under hypoxic hypothermia.

Snow finches (Onychostruthus taczanowskii, & Pyrgilauda ruficollis) and tree sparrow (Passer montanus) are closely related within the Old World Sparrows (Passeridae). Snow finches are endemic species on the Qinghai-Tibet Plateau (QTP) with a strict elevation distribution rang from 3500 to 5100m , whereas the tree sparrow is an introduced colonizer on the QTP with altitudinal distribution from sea level to 4400m (Figure 1 a) [8]. Snow finches and high-altitude tree sparrow separately have relative longer and 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 three objectives in this study. First, we found the 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 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.

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 (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 or myofibril (Figure 1b and Table 1). In contrast to the increased muscle-fiber size in highland tree sparrows, snow finches had thicker myofibril (Table 1), indicating that different modifications of muscle phenotype occur in intraspecies (high-altitude tree sparrow vs. low-altitude tree sparrow) and interspecies (two snow finches vs. low-altitude tree sparrow). Additionally, the ratio of wing length to body length significantly associated with muscle fiber area (R²=0.36, p= 0.0002, Figure S6a) and myofibril diameter (R²=0.51, p< 0.0001, Figure S6b).

The improvement in oxidative capacity of fight muscle 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 (Figure 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 (Figure 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²=0.44, p= 0.0004, Figure S6c) and myofibril size (R²=0.56, p< 0.0001, Figure S6d). 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 fight muscle was greater in snow finches than highland tree sparrows.

The first three principal components of the PCA on the basis of 12 phenotypic traits explained most of variations in pectoralis across inter- and intra-species (Figure S1). Within species, the PC1 (59.1%) represented mainly most of 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, whilst PC2 (18.4%) and PC3 (6.5%) explained muscle fiber density, capillary size as well as lipid droplet, respectively.

Gene expression differences and WGCNA

Transcriptomic analysis identified 966 differentially expressed genes (DEGs) in highland tree sparrow and 2457 DEGs shared in two snow finches comparing to lowland tree sparrow (Table S1). 407 DEGs, comprising 196 upregulated genes and 211 downregulated genes, appeared to all highland birds, enriched for “muscle structure development” (P < 0.001) and “oxygen transport” (P=0.0303) (Table S2).

Weighted gene correlation network analysis using WGCNA organized 966 DEGs within populations and 2457 shared DEGs in snow finches to four and five modules respectively (Figure 2a, b; Figure S2-S4). In population network (high- and low-altitude tree sparrows), blue module and brown module correlated positively with PC1 (Figure 2c). The DEGs in blue module were enriched into “muscle structure development” (P <0.001) and “muscle contraction” (P =0.0281) (Table S3). 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 S3). Also, yellow module and turquoise module associated negatively with PC1 (Figure 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 S3). As for species network, blue module and turquoise module associated positively with PC1 (Figure 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 “inositol phosphate metabolism” (P =0.001) (Table S3). And yellow module and brown module associated negatively with PC1 (Figure 2d). The DEGs in these modules were respectively enriched into “organelle” (P =0.009) and “membrane-bounded organelle” (P < 0.001) (Table S3). Only green module had a negative correlation with PC2 (Figure 2b) and genes in this module were related to “membrane part” (P =0.004) (Table S3). 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.

Discovery of candidate genes underlying phenotypic variations

The hub genes in co-expression networks are considered having the potential biological relevance and importance in driving phenotypic variations, so we identified hub genes within each module by connectivity and degree (top 30% and 15% of genes as hub genes in population and species network, respectively) (Table S4). As for population network, 9 hub genes in blue module involved in blood vessel formation (VEGFD, EGF, MMRN2), muscle cell differentiation (MYOZ3, PPAC3, SMRD3, FBX40), mitochondrial distribution (CLUH) and mitophagy (BNIP3) (Figure 2e; Figure 3a; Table S4). And eight genes were shared DEGs with snow finches, pointing to a common genetic basis of muscle phenotype. Of these genes, CLUH positively associated with DS within altitudinal populations (R²=0.87, p< 0.0001) but not species (R²=0.15, p= 0.1102) (Table S5). BNIP3 acted as a pro-apoptotic factor combining with autophagosome [9], and its expression connected negatively with the Vv(mt) within altitudinal populations (R²=0.41, p< 0.024) and species (R²=0.65, p< 0.0001) (Figure 3b, Table S5). 11 hub genes in turquoise module were involved in blood vessel development (PGS2, KPCA, NDST1, FLVC1, SEMA3C), muscle system process (MEF2C, BTG1, PRIC1), and mitochondrion localization (TRAK1, WIPI2, TRAK2) (Figure 2e; Figure 3a; Table S4). SEMA3C affected endothelial cell proliferation through increasing integrin activity similar to those induced by VEGF [10], and its expression correlated positively with capillary area within altitudinal populations (R²=0.64, p= 0.002) (Table S5, Figure 3b). MEF2C was related to striated muscle hypertrophy and obviously connected with fiber area within altitudinal populations (Figure 3b; Table S5). TRAK2 enhanced mitochondrial transport in axons of hippocampal neurons [11], and demonstrated a positive correlation with the DS within altitudinal populations (R²=0.67, p= 0.001) and species (R²=0.73, p< 0.0001) (Figure 3b; Table S5). Additionally, WIPI2 was suggested a mitophagy factor [9], and its expression showed a negative association with the Vv(mt) (Figure 3b; Table S5).

As for species network, two hub genes of the blue module were involved in vasculature development (MEF2C), and muscle structure development (MEF2C, ASB2) (Figure 2f; Figure 3c; Table S4). The turquoise module had 12 hub genes, 11 of these (expect for WIPI2) only in snow finches, probably indicating the regulatory difference of pectoralis phenotypes between snow finches and highland tree sparrow. The five genes in this 12 hub gene subset were involved in angiogenesis (EPAS1, VEGFR3, VEGFR4, PIK3C2A, KIAA1462), three were involved in muscle structure development (FLNB, MAML1, EPAS1), and six were involved in mitochondrion organization (LIG3, HIGD2A, WIPI2, EPAS1, MID51, HTT) (Figure 2f; Figure 3c; Table S4). 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 (Figure 2f; Table S5) [12-14]. 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²=0.28, P=0.0074) and CF (VEGFR4, R²=0.62, P<0.0001) across species (Figure 3c, Table S5). MID51 involving in mitochondrial fission was a hub gene with a higher degree and its expression positively associated with the Vv(mt) (R²=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 (Figure 2f; Figure 3a; Table S5).

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 (Figure 4a, b; Figure S5a), indicating improvement in glucose uptake and utilization. Although plasma insulin concentration was similar to those of lowland birds (Figure 4a), highland birds had some DEGs involving in regulating insulin sensitivity and insulin signaling pathway (AR/ERα, PIK3C2A, PIK3C2B, SLC2A12 (Glut12) and PTP1B) (Figure 4c). Additionally, expression and activities of some glycolytic enzymes (hexokinase, pyruvate kinase, and lactate dehydrogenase) were increased in highland birds (Figure 4a, c). Consistently, glycolytic intermediates increased in all highland birds, with increased glucose and fructose-6-phosphate (Figure 4b). Activity of citrate synthase (CS) was higher in highland birds (Figure 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 (Figure 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 two snow finches (Figure 5a; Figure S5b). The ratio of short-chain fatty acid (SFA), however, was great in snow finches instead of a low ratio of desaturated fatty acid (DFA) (Figure 5a). Meanwhile, the ratio of long-chain fatty acid (LFA) also significantly decreased in snow finches (Figure 5a), indicating DFA and LFA might be important fuel substrates. Additionally, PPARδ had a low expression level in snow finches, whilst DGAT2 activating lipid depot formation was up-regulated and its expression was positively associated with volume density of lipid droplet (R²=0.49, P= 0.0002) (Figure 5b; Table S5). 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 (Figure 5a, b, c). Together, snow finches likely decrease PPARδ expression to promote intramuscular lipid biosynthesis and increase mitochondrial metabolism for fatty acid oxidation in a long evolutionary time to highland. In contrast, hydroxybutyrate (HB) decreased significantly in highland tree sparrow (Figure 5a), indicating that ketone instead of fatty acid was a major anaplerotic substrate for energetic supply in a short evolutionary time to hypoxia.

By integrating a comprehensive array of analytical techniques, we illuminate muscle phenotypic variation and its underpinning regulatory as well as metabolic mechanism across altitudinal songbirds with different evolutionary time. Despite considerable increases in muscle oxidative capacity, capillarity and mitochondrial abundance when facing hypothermia and hypoxia, different regulatory mechanisms might drive their occurrence in snow finches and high-altitude tree sparrow. A metabolic feature shared by 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 for different evolutionary time under hypothermia.

Pectoralis variations to highland survival

Great body mass and muscle weight in highland animals are concordance with Bergmann’s rule[15, 16]. Avian pectoralis is used to support the high metabolic costs through aerobic metabolism for flight or thermogenesis [4]. High-altitude birds, however, suffer from severely limited oxygen 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 [2]. Either size modification of muscle fiber or myofibril in highland birds plays an important role on heat production through shivering thermogenesis for maintenance of a constant core body temperature [17]. Various modifications of muscle phenotype across altitudinal birds are probably associated with the flight performance in hypoxic environment for an increased size of myofibril which would produce more contractile force to prevent flightless, probably accounting for why highland tree sparrows perform a weak flight capacity [18].

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 [19]. Capillary size increase may reduce capillary flow resistance directly and maintain capillary rheology in hypoxic environment [20]. Mitochondrion is the main place of aerobic oxidation and takes a majority of ATP production in eukaryote [21]. An increase in DS in highland birds may enhance ability to produce proton-motive force as reported in previous studies [22]. Overall, capacities of oxidation and oxygen delivery in pectoralis are prompted through differentially phenotypic variations, probably resorting from different adaptation time to hypoxic hypothermia and greater Vv (mt) in snow finches probably contributes to fat metabolism to promote myofibril contraction and maintain thermogenesis or flight.

Potential transcriptomic basis of pectoralis variations

Given that the differences of 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 regulation of muscle development via MRF4-MEF2 axis, which is thought to a novel pathway regulating muscle fiber size and muscle mass in adult muscle tissue [23]. MRF4 and MEF2C are known to play a role in myogenesis [24]. MRF4 knockdown in skeletal muscle of adult rodent stimulates MEF2 transcriptional activity and then causes muscle hypertrophy [23]. All highland birds promote the expression of MEF2C and its target genes (i.e. MYH1, ACTA; Table S1; Figure 2c; Figure 3a), but only highland tree sparrow show the greater fiber size (Table 1). Instead, snow finches promote a muscle hyperplasia with greater myofibril. The overexpression of hub gene EPAS1 likely promotes HDAC3 expression which can deacetylate and repress MEF2C, and thus causes muscle hyperplasia as previously reported in mice [25]. 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 is also a regulative process of postnatal but not embryonic myogenesis [23].

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 [10]. Additionally, snow finches also have an increase in capillary density along with the overexpression of EPAS1 and its downstream genes (i.e. VEGFR3/4, DLL4, HEY2, TIE2; Figure 3a). Mounting evidence suggests that EPAS1 has an important role in angiogenesis though most of studies are focused on HIF1α [12-14]. 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 [12-14]. 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 [26]. Nevertheless, the up-regulated expression of TIE2 and ANGPTL1 stimulates basement membrane deposition and pericytes detachment, thereby mediating vessel maturation [27].

Autophagy is a lysosome-mediated degradation process for cytoplasmic components, and also plays an important role in mitochondrial homeostasis [28]. Highland birds decrease the expression of genes involved in mitophagy, including WIPI2, LC3, and BNIP3 (Figure 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 [11]. However, gene expression for mitochondrial fission and fusion is various. Snow finches and highland tree sparrow induce mitochondrial fission via MID51 and MID49, respectively [9]. Additionally, MIGA1 and MIGA2 also enhance mitochondrial fusion only in snow finches [29]. Besides, KLF4/ERR/PGC1a transcriptional complex regulates a broad spectrum of genes involved in mitochondrial biogenesis, dynamics, and metabolism [30]. All three 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 prevailing views considers energetic needs for thermogenesis in highland coming from glucose oxidation for which requires less O₂ than fatty acid oxidation [31]. Glucose uptake from the bloodstream into cells is mediated by a family of facilitative glucose transporters (Gluts) [6]. It has been long suggested a low sensitivity to insulin due to the absence of Glut 4 in birds [7]. 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 [32]. Highland birds exhibit a similar concentration of plasma insulin and an accelerated glucose uptake in pectoralis major along with high expression of Glut12, indicating an improvement of insulin sensitivity. Previous studies have suggested that AR/ ERα, PIK3C2A, PIK3C2B and PTP1B mediate insulin sensitivity through activating insulin signaling pathway [33]. It is the first time to find an enhanced glucose metabolism via activating insulin signaling pathway by gene expression in highland vertebrates, differing from through activity change of PTP1B in insects [34]. Additionally, increase in insulin sensitivity improves glucose oxidation[35], which mainly exhibits in increased capacities of glycolysis and TCA cycle [36].

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 [31], thereby it is an important energy strategy taken by small animal when cold exposure [37]. The enhanced capacities of lipid oxidation probably occurs in snow finches facing prolonged cold threaten 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 [36, 37].

Determining the molecular and functional mechanisms underlying phenotypic variation is a fundamental goal of evolutionary biology and ecology. Pectoralis is a fundamental tissue for locomotion and thermogenesis and its variation is very vital for bird survival in highland. However, the regulatory and metabolic basis of pectoralis variation remains largely unknown for highland birds, in particular resident species with different evolutionary history. In this study, we find that highland birds have an extensive increase in capillarity, mitochondrial abundance and oxidative capacity of muscle ﬁber. RNA-seq analyses reveal that 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 drives their occurrence in part due to own various evolutionary histories. Moreover, a shared metabolic feature among 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. 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.

Ethics statement

This study was carried out in strict accordance with the recommendations of the Regulations for the Administration of Affairs Concerning Experimental Animals (Ministry of Science and Technology, China, revised in June 2004). All procedures performed on birds were approved by the Institute of Zoology Animal Care Committee. All birds were euthanized through cardiac compression, and all efforts were made to minimize suffering of animals.

Sampling

Wild birds were caught with mist-net in summer of 2016, highland bird populations from Qinghai Province about 3200m (10 tree sparrows, Passer montanus) and 3900m (10 rufous-necked snow finch, Pyrgilauda ruficollis and 10 white-rumped snow finch, Onychostruthus taczanowskii) and lowland population (n= 8) of the tree sparrow from Yanqi Lake of Beijing at 80m (many details were represented in table S6). All birds were euthanized in 3-5 seconds through cardiac compression. Pectoral major muscle was dissected immediately following euthanasia, flash-frozen in liquid nitrogen and stored at -80°C.

Muscle histology and transmission electron microscopy

Oxidative muscle type and capillarity were performed as previously described [38]. Brieﬂy, 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 N2). Muscle sections (10 mm) were obtained in a Cryostat Microtome (Leica CM900, Germany) at -20℃. Oxidative muscle ﬁbers were identiﬁed by succinate dehydrogenase activity, by staining in assay buffer (concentrations in mM: 0.6 nitroblue tetrazolium, 2.0 KH₂PO₄, 15.4 Na₂HPO₄, 16.7 sodium succinate) for 1 h at room temperature. Using alkaline phosphatase activity identified muscle capillaries, also by staining for 1 h at room temperature (assay buffer concentrations in mM: 1.0 nitroblue tetrazolium, 0.5 5-bromo-4-chloro-3-in-doxyl phosphate, 28 NaBO2,7 MgSO₄; pH 9.4). The sections were imaged using light microscopy and sufﬁcient images (eight or more) were analysed for each sample using image J software. Biochemicals were obtained from Sigma-Aldrich (Shanghai, China).

The fight muscle was removed from an intermediate depth and was then ﬁxed at 4°C for 24–48 h in 2% glutaraldehyde in 0.1M PBS buffer at pH 7.4. Small muscle blocks (2mm x 2mm) were prepared and post-ﬁxed 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 sectionswere 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 [39]. Grid size of 90 nm was used to estimate mitochondria and lipid droplet volume at a square size of 4460*4460nm.

Plasma glucose, lactate, insulin, tissue glycogen measurement, and enzyme activity assays

To quantify circulating glucose levels, fasting plasma was obtained from eight birds from each species and was measured using an Accu-Check 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.

The pectoralis was homogenized in 10 volumes of ice-cold homogenization buffer (in mM: 50 hepes, 5 EDTA, 0.2 dithiothreitol (DTT), and 0.1% Triton-X-100; pH 7.4). After homogenates were centrifuged at 1000g at 4℃, the supernatants were collected and determined for protein contents using the G-250 method. The maximal activities of 5 enzymes (hexokinase, pyruvate kinase, lactate dehydrogenase, 3-hydroxyacyl-CoA dehydrogenase, and citrate synthase) in pectoralis were measured at sparrow rectal temperature (40℃) in 100mmol l^-1 KH₂PO₄ (pH 7.4) with a 96-well plate ELIASA. All assays were optimized to assure substrates and enzymes were not limited and carried out in triplicate with the following reaction conditions (in mM, unless stated otherwise).

For carbohydrate metabolism: (i) hexokinase (HK): 10 mmol l⁻¹ glucose, 3 mmol l⁻¹ ATP, 10 mmol l⁻¹ MgCl₂, 1.5 mmol l⁻¹ NADP⁺, 1 unit of glucose-6-phosphate dehydrogenase; (ii) pyruvate kinase (PK): 10mmoll⁻¹ phosphoenol pyruvate, 2.5 mmol l⁻¹ ADP, 10 mmol l⁻¹ MgCl_2, 0.15 mmol l⁻¹ NADH, 1 unit of lactate dehydrogenase; (iii) lactate dehydrogenase (LDH): 5 mmol l⁻¹ pyruvate, 0.15 mmol l⁻¹ NADH. For fatty acid metabolism: (i) 3-hydroxyacyl-CoA dehydrogenase (HOAD): 0.1 mmol l⁻¹ acetoacetyl CoA, 0.5 mmol l⁻¹ NADH, 0.2 DTT. For the tricarboxylic acid cycle: (i) citrate synthase (CS): 0.5 mmol l⁻¹ oxaloacetate, 0.15 mmol l⁻¹ acetyl-coA, 0.15 mmol l⁻¹ 5,5′- dithiobis-2-nitrobenzoic acid, 0.1% Triton X-100, in 100 Tris (pH 8.0). Background rates were calculated in control reactions lacking specific substrates. Extinction coefficients were 6.22 (mmol l⁻¹)⁻¹cm⁻¹for NADH (340nm) and 13.6 (mmol l⁻¹)⁻¹cm⁻¹for DTNB (412nm) to measure enzyme activities.

RNA isolation, library preparation, sequencing, quantification, and normalization

Total RNA was extracted from each 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. Briefly, mRNA was purified from total RNA using poly-T oligo-attached magnetic beads. Fragmentation was carried out using divalent cations under elevated temperature in NEB Next First Strand Synthesis Reaction Buffer(5X). First strand cDNA was synthesized using random hexamer primer and M-MuLV. Reverse Transcriptase(RNase H-). Second strand cDNA synthesis was subsequently performed using DNA Polymerase I and RNase H. Remaining overhangs were converted into blunt ends via exonuclease/polymerase activities. After adenylation of 3’ ends of DNA fragments, NEBNext Adaptor with hairpin loop structure were ligated to prepare for hybridization. In order to select cDNA fragments of preferentially 250~300 bp in length, the library fragments were purified with AMPure XP system (Beckman Coulter, Beverly, USA). Then 3 µl USER Enzyme (NEB, USA) was used with size-selected, adaptor-ligated cDNA at 37°C for 15 min followed by 5 min at 95 °C before PCR. Then PCR was performed with Phusion High-Fidelity DNA polymerase, Universal PCR primers and Index (X) Primer. At last, PCR products were purified (AMPure XP system) and library quality was assessed on the Agilent Bioanalyzer 2100 system.

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 [40] was used to filter reads containing adapter, reads containing ploy-N and low quality reads based on read quality checked with FASTQC [41]. The parameters used were as follows: sliding window = 4-bp; Phred33 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 (unpublished data) using STAR with default parameters [42]. 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 three species were obtained by intersecting the lists of above two orthologs. After the reads were mapped to the reference genomes [43], expression quantifications of genes were performed using RSEM [44]. Only genes with count number greater than 1 in at least four samples were included in differential gene expression analysis. Expression levels for genes with one-to-one orthologues in all three bird species (n = 12951) were normalized with a RLE (relative log expression) method across muscle samples [45].

Differential expression analysis and weighted gene coexpression network analysis (WGCNA)

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 fltering was enabled in a R/Bioconductor package Deseq2 (R version 3.5.1) with a false discovery rate (FDR) < 0.05 based on Benjamini-Hochberg method to control the False Discovery Rate in multiple tests context in identifying significantly differentially expressed genes [46]. The cut off values for log-2-fold change were set at 0.59 and -0.59.

And 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, log2-transformed counts to analyze the DEGs with the blockwiseModules function in WGCNA [47] 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 [48]. 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 using a principal component analysis (PCA) of gene expression with the blockwiseModules 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 bub 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) [49]. We then implemented gene ontology categories within modules positively and negatively correlated with muscle traits using G:PROFILER (FDR < 0.05) [50]. According to the results of functional GO enrichment analysis, intramodular hub genes were identified as candidate genes driving potentially muscle variation.

Gas chromatography-mass spectrometry analysis

Pectoralis (10 lowland tree sparrows, 10 highland tree sparrows, 10 rufous-necked snow finches, 10 white-rumped snow finches) and plasma (7 lowland tree sparrows, 7 highland tree sparrows, 6 rufous-necked snow finches, 7 white-rumped snow finches) metabolites were extracted by methanol/chloroform protocol, as described previously [51]. 1 ml of chloroform/methanol/ distilled water mixture (chloroform: methanol: distilled water = 1:4:1 mixture) was added to Eppendorf tubes containing ∼50 mg pectoralis powdered in lipid nitrogen. 450 μL of methanol/distilled water mixture (8:1 mixture) was added to Eppendorf tubes containing ∼50 uL plasma. 20 μg, and 10 μg of heptadecanoic acid as well as decanoic acid as double internal standard were added into muscle and plasma tube, respectively. After vortex for 10 s, samples were kept on ice for 15 min and then in a sonication bath for 15 min. After centrifugation at 12,000 rpm for 15 min at 4 °C, 200 μL of supernatant was transferred to a 2 mL auto-sampler vial and dried in a vacuum oven. Dried samples were derivatized using methoxyamine hydrochloride solution (50 μL of 15mg·mL−1 methoxyamine in pyridine). Themixture was kept for 16 h at room temperature for methoxymation, and then 50 μL of TMSFA containing 1% TMCS was added for trimethylsilylated. After 1 h trimethysilylation, 30 μL of hexane was added and then transferred to an insert in 2 mL autosampler vial for GC–MS analysis. A 2 μL of derivatized sample was injected by an Agilent 4683B GC auto-sampler (Agilent Technologies, Atlanta, USA) into an Agilent 6890N gas chromatograph with 5973 mass spectrometry at 250 °C without split. The detail program was set as previous study [50]. GC-MS chromatograms were processed using GAVIN based on Matlab [52]. Deconvolution was achieved by AMDIS and the results were imported into GAVIN for retention time correction and area integration across the data set. Identification of metabolites from GC-MS analysis was supported by comparison with a standard mix (Supelco 37 Component FAME Mix; Sigma-Aldrich; six monosaccharide quantitative standards; Ludger).

Phylogenetic analyses

We used MrBayes 3.2 [53] for the phylogenetic tree construction with the whole mitochondrial genome using best-fit models based on BIC model selection criteria (Figure 1 a). Sequences were downloaded from GenBank (Accession: NC_024821.1, NC_022815.1 and NC_025914.1 for Passer montanus, Pyrgilauda ruficollis and Onychostruthus taczanowskii, respectively).

Statistical analyses

The two-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 were considered significant.

GC-MS: Gas chromatography-mass spectrometry

DEGs: Differentially expressed genes

WGCNA: Weighted gene coexpression network analysis

GO: Gene Ontology;

KEGG: The Kyoto Encyclopedia of Genes and Genomes

Acknowledgements

Many thanks to Dehua Wang, Zhiyun Jia, Graham R. Scott, Shane G. DuBay, Qingsheng Chi, Boyang Ding and Ding Ding for discussion with laboratory techniques. Thanks also to Zuohua Yin, Dan Zhu, Yinzi Ma, and Pengyan Xia for fieldwork, laboratory space, and transmission electron microscopy.

Funding

This work was funded by National Natural Science Foundation of China (No. 31630069 to F.L.), the Strategic Priority Research Program of the Chinese Academy of Sciences (Grant Nos. XDA19050202 and XDB13020300 to F.L.) and the Second Tibetan Plateau Scientific Expedition and Research (STEP) program (Grant No. 2019QZKK0501 to F.L.).

Availability of data and materials

Raw sequencing data have been deposited at the National Center for Biotechnology Information short read archive under bioproject PRJNA523449.

Authors’ contributions

FML and YX designed and performed the experiments and wrote the manuscript. YH, YLC, LQF, GS and DML helped bioinformatics analysis and fieldworks for sample collection. YHQ gave many comments for the manuscript, and all authors commented and approved.

Ethics approval and consent to participate

All procedures performed on animal were approved by the institute of zoology animal care committee.

Competing interests

The authors of this manuscript have declared no competing interests.

Consent for publication

Not applicable

Author details

¹Key Laboratory of Zoological Systematics and Evolution, Institute of Zoology, Chinese Academy of Sciences, Beijing, 100101, China. ²University of Chinese Academy of Sciences, Beijing, 100049, China. ³Key Laboratory of Forest Ecology in Tibet Plateau of Ministry of Education, Institute of Tibet Plateau Ecology, Tibet Agriculture & Animal Husbandry University, Nyingchi 860000, China. ⁴Key Laboratory of Animal Physiology, Biochemistry and Molecular Biology of Hebei Province, College of Life Sciences, Hebei Normal University, Shijiazhuang 050024, P. R. China. ⁵Center for Excellence in Animal Evolution and Genetics, Chinese Academy of Sciences, Kunming, 650223, China

Scott GR, Elogio TS, Lui MA, Storz JF, Cheviron ZA: Adaptive modifications of muscle phenotype in high-altitude deer mice are associated with evolved changes in gene regulation. Mol Biol Evol 2015, 32(8):1962-1976.
Scott GR, Egginton S, Richards JG, Milsom WK: Evolution of muscle phenotype for extreme high altitude flight in the bar-headed goose. P R Soc B 2009, 276(1673):3645-3653.
Schiaffino S, Reggiani C: Fiber types in mammalian skeletal muscles. Physiol Rev 2011, 91(4):1447-1531.
Marquez J, Sweazea KL, Braun EJ: Skeletal muscle fiber composition of the English sparrow (Passer domesticus). Comp Biochem Phys B 2006, 143(1):126-131.
Qu YH, Zhao HW, Han NJ, Zhou GY, Song G, Gao B, Tian SL, Zhang JB, Zhang RY, Meng XH et al: Ground tit genome reveals avian adaptation to living at high altitudes in the Tibetan plateau. Nat Commun 2013, 4.
Shepherd PR, Kahn BB: Glucose transporters and insulin action—implications for insulin resistance and diabetes mellitus. New Engl J Med 1999, 341(4):248-257.
Polakof S, Mommsen TP, Soengas JL: Glucosensing and glucose homeostasis: from fish to mammals. Comp Biochem Phys B 2011, 160(4):123-149.
del Hoyo J, Elliot A, Sartagal J, Christie D, De Juana E: Handbook of the birds of the world alive. Barcelona: Lynx Edicions. 2014. In.; 2018.
Ni HM, Williams JA, Ding WX: Mitochondrial dynamics and mitochondrial quality control. Redox Biol 2015, 4:6-13.
Banu N, Teichman J, Dunlap-Brown M, Villegas G, Tufro A: Semaphorin 3C regulates endothelial cell function by increasing integrin activity. Faseb J 2006, 20(12):2150-2152.
MacAskill AF, Kittler JT: Control of mitochondrial transport and localization in neurons. Trends Cell Biol 2010, 20(2):102-112.
Rankin EB, Rha J, Unger TL, Wu CH, Shutt HP, Johnson RS, Simon MC, Keith B, Haase VH: Hypoxia-inducible factor-2 regulates vascular tumorigenesis in mice. Oncogene 2008, 27(40):5354-5358.
Skuli N, Liu LP, Runge A, Wang T, Yuan LJ, Patel S, Iruela-Arispe L, Simon MC, Keith B: Endothelial deletion of hypoxia-inducible factor-2 alpha (HIF-2 alpha) alters vascular function and tumor angiogenesis. Blood 2009, 114(2):469-477.
Park YS, Kim G, Jin YM, Lee JY, Shin JW, Jo I: Expression of angiopoietin-1 in hypoxic pericytes: Regulation by hypoxia-inducible factor-2 alpha and participation in endothelial cell migration and tube formation. Biochem Bioph Res Co 2016, 469(2):263-269.
Fan L, Cai T, Xiong Y, Song G, Lei F: Bergmann’s rule and Allen’s rule in two passerine birds in China. Avian Research 2019, 10(1):34.
Sun Y, Li M, Song G, Lei F, Li D, Wu Y: The role of climate factors in geographic variation in body mass and wing length in a passerine bird. Avian Research 2017, 8(1):1.
Chaillou T: Skeletal muscle fiber type in hypoxia: adaptation to high-altitude exposure and under conditions of pathological hypoxia. Front Physiol 2018, 9.
Sun YF, Ren ZP, Wu YF, Lei FM, Dudley R, Li DM: Flying high: limits to flight performance by sparrows on the Qinghai-Tibet Plateau. J Exp Biol 2016, 219(22):3642-3648.
Scott GR, Milsom WK: Control of breathing and adaptation to high altitude in the bar-headed goose. Am J Physiol-Reg I 2007, 293(1):R379-R391.
Hauck EF, Apostel S, Hoffmann JF, Heimann A, Kempski O: Capillary flow and diameter changes during reperfusion after global cerebral ischemia studied by intravital video microscopy. J Cerebr Blood F Met 2004, 24(4):383-391.
Rich P: Chemiosmotic coupling: the cost of living. Nature 2003, 421(6923):583.
Mahalingam S, McClelland GB, Scott GR: Evolved changes in the intracellular distribution and physiology of muscle mitochondria in high-altitude native deer mice. J Physiol-London 2017, 595(14):4785-4801.
Moretti I, Ciciliot S, Dyar KA, Abraham R, Murgia M, Agatea L, Akimoto T, Bicciato S, Forcato M, Pierre P et al: MRF4 negatively regulates adult skeletal muscle growth by repressing MEF2 activity. Nat Commun 2016, 7.
Hernandez-Hernandez M, Garcia-Gonzalez EG, Brun CE, Rudnicki MA: The myogenic regulatory factors, determinants of muscle development, cell identity and regeneration. Semin Cell Dev Biol 2017, 72:10-18.
Wu MZ, Tsai YP, Yang MH, Huang CH, Chang SY, Chang CC, Teng SC, Wu KJ: Interplay between HDAC3 and WDR5 is essential for hypoxia-induced epithelial-mesenchymal transition. Mol Cell 2011, 43(5):811-822.
Hellstrom M, Phng LK, Hofmann JJ, Wallgard E, Coultas L, Lindblom P, Alva J, Nilsson AK, Karlsson L, Gaiano N et al: Dll4 signalling through Notch1 regulates formation of tip cells during angiogenesis. Nature 2007, 445(7129):776-780.
Augustin HG, Koh GY, Thurston G, Alitalo K: Control of vascular morphogenesis and homeostasis through the angiopoietin-Tie system. Nat Rev Mol Cell Bio 2009, 10(3):165-177.
Lee J, Giordano S, Zhang JH: Autophagy, mitochondria and oxidative stress: cross-talk and redox signalling. Biochem J 2012, 441:523-540.
Zhang YP, Liu XM, Bai J, Tian XJ, Zhao XC, Liu W, Duan XY, Shang WN, Fan HY, Tong C: Mitoguardin regulates mitochondrial fusion through MitoPLD and is required for neuronal homeostasis. Mol Cell 2016, 61(1):111-124.
Dorn GW, Vega RB, Kelly DP: Mitochondrial biogenesis and dynamics in the developing and diseased heart. Gene Dev 2015, 29(19):1981-1991.
Schippers MP, Ramirez O, Arana M, Pinedo-Bernal P, McClelland GB: Increase in carbohydrate utilization in high-altitude andean mice. Curr Biol 2012, 22(24):2350-2354.
Coudert E, Pascal G, Dupont J, Simon J, Cailleau-Audouin E, Crochet S, Duclos MJ, Tesseraud S, Metayer-Coustard S: Phylogenesis and biological characterization of a new glucose transporter in the chicken (Gallus gallus), GLUT12. Plos One 2015, 10(10).
Yu IC, Lin HY, Liu NC, Sparks JD, Yeh SY, Fang LY, Chen LM, Chang CS: Neuronal androgen receptor regulates insulin sensitivity via suppression of hypothalamic nf-kappa b-mediated ptp1b expression. Diabetes 2013, 62(2):411-423.
Ding D, Liu GJ, Hou L, Gui WY, Chen B, Kang L: Genetic variation in PTPN1 contributes to metabolic adaptation to high-altitude hypoxia in Tibetan migratory locusts. Nat Commun 2018, 9.
Yan Lj: Redox imbalance stress in diabetes mellitus: Role of the polyol pathway. Animal models and experimental medicine 2018, 1(1):7-13.
Horscroft JA, Kotwica AO, Laner V, West JA, Hennis PJ, Levett DZ, Howard DJ, Fernandez BO, Burgess SL, Ament Z: Metabolic basis to Sherpa altitude adaptation. Proceedings of the National Academy of Sciences 2017, 114(24):6382-6387.
Cheviron ZA, Bachman GC, Connaty AD, McClelland GB, Storz JF: Regulatory changes contribute to the adaptive enhancement of thermogenic capacity in high-altitude deer mice. P Natl Acad Sci USA 2012, 109(22):8635-8640.
Scott GR, Johnston IA: Temperature during embryonic development has persistent effects on thermal acclimation capacity in zebrafish. P Natl Acad Sci USA 2012, 109(35):14247-14252.
Weibel E: Profile size and particle size Weibel ER eds.. Stereological methods: practical methods for biological morphometry, Vol. 1: 51-60. In.: Academic Press, Inc. New York; 1979.
Bolger AM, Lohse M, Usadel B: Trimmomatic: a flexible trimmer for Illumina sequence data. Bioinformatics 2014, 30(15):2114-2120.
Andrews S: FastQC: a quality control tool for high throughput sequence data. 2010. In.; 2016.
Dobin A, Davis CA, Schlesinger F, Drenkow J, Zaleski C, Jha S, Batut P, Chaisson M, Gingeras TR: STAR: ultrafast universal RNA-seq aligner. Bioinformatics 2013, 29(1):15-21.
Qu Y, Chen C, Xiong Y, She H, Zhang Y, Cheng Y, Dubay SG, Li D, Ericson PGP, Hao Y: Rapid phenotypic evolution with shallow genomic differentiation during early stages of high elevation adaptation in Eurasian Tree Sparrows. Natl Sci Rev 2020, 7(1):113-127.
Li B, Dewey CN: RSEM: accurate transcript quantification from RNA-Seq data with or without a reference genome. BMC Bioinformatics 2011, 12(1):323.
Dillies M, Rau A, Aubert J, Hennequetantier C, Jeanmougin M, Servant N, Keime C, Marot G, Castel D, Estelle J: A comprehensive evaluation of normalization methods for Illumina high-throughput RNA sequencing data analysis. Briefings in Bioinformatics 2013, 14(6):671-683.
Love MI, Huber W, Anders S: Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol 2014, 15(12).
Langfelder P, Zhang B, Horvath S: Defining clusters from a hierarchical cluster tree: the Dynamic Tree Cut package for R. Bioinformatics 2008, 24(5):719-720.
Zhang B, Horvath S: A general framework for weighted gene co-expression network analysis. Stat Appl Genet Mo B 2005, 4.
Han JDJ, Bertin N, Hao T, Goldberg DS, Berriz GF, Zhang LV, Dupuy D, Walhout AJM, Cusick ME, Roth FP et al: Evidence for dynamically organized modularity in the yeast protein-protein interaction network. Nature 2004, 430(6995):88-93.
Reimand J, Kull M, Peterson H, Hansen J, Vilo J: g : Profiler - a web-based toolset for functional profiling of gene lists from large-scale experiments. Nucleic Acids Res 2007, 35:W193-W200.
Shi YL, Chi QS, Liu W, Fu HP, Wang DH: Environmental metabolomics reveal geographic variation in aerobic metabolism and metabolic substrates in Mongolian gerbils (Meriones unguiculatus). Comp Biochem Phys D 2015, 14:42-52.
Behrends V, Tredwell GD, Bundy JG: A software complement to AMDIS for processing GC-MS metabolomic data. Anal Biochem 2011, 415(2):206-208.
Ronquist F, Huelsenbeck JP. MrBayes 3: Bayesian phylogenetic inference under mixed models. Bioinformatics 2003; 19:1572-1574.

Table 1 Masses and phenotypes of pectoralis in sparrows

	*Pa. mo* (L)	*Pa. mo* (H)	*Pa. ru*	*On. ta*
Organ mass (mg g^-1 body mass)^a
PMM (g)	1.21±0.04	1.91±0.07*	1.73±0.11*	1.90±0.10*
RM_PM (%)	7.50±0.18	8.38±0.30*	6.66±0.23†	5.45±0.48†
Fiber type and capillary density in pectoralis ^a
FD (mm^-2)	1723.20±44.79	1373.94±41.74†	1887.86±28.06*	1768.51±33.67
FS (um²)	583.12±15.46	733.09±19.82*	530.74±7.77†	567.27±10.66
MD(nm)	897.66±41.30	921.90±35.49	1134.71±30.98*	1122.51±40.14*
CD(mm^-2)	2869.12±69.68	2506.09±65.23†	3731.58±93.67*	3409.45±67.62*
CA(um^-2)	8.65 ±0.69	12.42±1.50*	13.26±1.50*	13.14±1.86*
CF	1.67±0.04	1.84±0.01*	1.98±0.03*	1.93±0.01*
Mitochondrial volume densities and numerical densities ^b
DS	0.314±0.024	0.481±0.021*	0.536±0.007*	0.518±0.013*
Vv (mt)	0.078 ±0.004	0.092 ±0.003*	0.114 ±0.004*	0.118 ±0.005*
Vv (ssm)	0.047 ±0.002	0.057 ±0.001*	0.069 ±0.002*	0.066 ±0.002*
Vv (imm)	0.032 ±0.003	0.034 ±0.001	0.046 ±0.003*	0.051 ±0.004*
Vv(sim)	0.015 ±0.002	0.023 ±0.001*	0.023 ±0.003*	0.015 ±0.003
Vv (LD)	0.0041±0.0006	0.0060±0.0014	0.0117±0.0017*	0.0133±0.0026*

NOTE.—Data are means ± standarad error ; * significant increase; † significant decrease. PMM, pectoralis major mass; RM_PM, ratio of pectoralis major mass to body mass; FD, fiber density; FS, fiber area; MD, myofibril diameter; CD, capillary density; CA, capillary area; CF, the number of capillaries per fiber; DS, proportion of subsarcolemmal mitochondrion; Vv (mt), total mitochondrial volume density; Vv (ssm), volume density of subsarcolemmal mitochondrion; Vv (imm), volume density of intermyofibrillar mitochondrion; Vv(sim), the difference between Vv (ssm) and Vv (imm); Vv (LD): volume density of lipid droplet. Pa. mo (L) (lowland tree sparrow), Pa. mo (H) (highland tree sparrow), Pa. ru (rufous-necked snow finch), On. ta (white-rumped snow finch).

a n=8, 10, 10 &10 for Pa. mo (L), Pa. mo (H), Pa. ru, On. ta.

b n=6 for each species.

FigureS1.pdf
Additional file 1: Figure S1. Rotated principal component analysis of pectoralis major. Variance between intra-species represented 59.1%, 13.1%, and 8% for PC1, PC2, and PC3, respectively. Variance between inter-species represented 58.7%, 18.4%, and 6.5% for PC1, PC2, and PC3, respectively. FD, fiber density; FS, fiber area; MD, myofibril diameter; CD, capillary density; CA, capillary area; CF, the number of capillaries per fiber; DS, proportion of subsarcolemmal mitochondrion; Vv .mt., total mitochondrial volume density; Vv .ssm., volume density of subsarcolemmal mitochondrion; Vv .imm., volume density of intermyofibrillar mitochondrion; Vv.sim., the difference between Vv .ssm. and Vv .imm.; Vv .LD.: volume density of lipid droplet.
TableS1.xlsx
Additional file 2: Table S1. Differentially expressed genes in altitudinal population and species pairwise analysis.
TableS2.docx
Additional file 3: Table S2. Gene ontology terms, Human Phenotype Ontology, and KEEG pathways enriched in shared differentially expressed genes.
FigureS2.pdf
Additional file 4: Figure S2. 18 and 22 of soft power for scalefree network construction of intraspecies (a) and interspecies (b) respectively. Correlation between co-expression modules and muscle traits across inter- (c) as well as intra-species (d).
FigureS3.pdf
Additional file 5: Figure S3. Modules correspond to muscle samples within species (a) and between species (b). Heat maps describing expression levels for all DEGs in muscle for each module, Barplots of the values of the module eigengene (i.e., the ﬁrst principal component) derived from singular value decomposition are displayed for each module.
FigureS4.pdf
Additional file 6: Figure S4. Network heatmaps of the 966 and 2457 genes differentially expressed between intraspecies (a) and interspecies (b) respectively. Pearson correlations between expression profiles of all pairs of genes were transformed into network connection strengths. Hub genes are highlighted in black squares.
TableS3.docx
Additional file 7: Table S3. Gene ontology terms, Human Phenotype Ontology, and KEEG pathways enriched in muscle functional module of differentially expressed genes.
TableS4.xlsx
Additional file 8: Table S4. Hub differentially expressed genes in modules associated with muscle phenotypes.
TableS5.docx
Additional file 9: Table S5. Candidate DEGs connected with phenotypic variations.
FigureS5.pdf
Additional file 10: Figure S5. Relative levels of intermediates and fatty acids in fasting plasma. * significant increase; ‡ significant decrease.
FigureS6.pdf
Additional file 11: 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).
TableS6.docx
Additional file 12: Table S6. Specimen information used for phenomic, transcriptomic and metabolomic analyses.

Download PDF

Journal Publication

published 21 Dec, 2021

Read the published version in Integrative Zoology →

Version 2

posted

You are reading this latest preprint version

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

Status:

Journal Publication

Version 2

Abstract

Figures

Background

Results

Discussion

Conclusions

Methods

Abbreviations

Declarations

References

Tables

Supplementary Files

Status:

Journal Publication

Version 2