Pectoral muscle transcriptome analyses reveal high-altitude adaptations in Tibetan chickens

The largest muscles in fowl are the pectorals, which provide the power required for birds to ﬂy. Tibetan chickens show speciﬁc adaptations to high-altitude conditions, but changes in the muscle transcriptome associated with these adaptations have not been characterized yet. Therefore, in this study, we used next-generation sequencing technologies to generate eight libraries of mRNA sequences for four Tibetan chickens and four Beijing fatty chickens. A comprehensive transcriptome analysis was performed. In the eight samples, 12333 annotated protein-coding genes were expressed. Among these, 48 differentially expressed genes were found; all of which were upregulated in Tibetan chickens. These differentially expressed genes were mainly involved in kidney morphogenesis, which indicates that hypoxia has an important effect on renal tubule development. Only nine genes were involved in Kyoto Encyclopedia of Genes and Genomes pathways, such as the endocytosis pathway, the MAPK signaling pathway, the calcium signaling pathway and the TGF-beta signaling pathway. The differentially expressed genes identiﬁed in this study will be used to facilitate future research into the Tibetan chicken.


Introduction
The molecular mechanisms underlying hypoxic adaptations in high-altitude animals have long garnered attention from biological and medical researchers . Adaptation to hypoxia in animals is a complex trait that involves multiple genes and multi-channel control. The genetic mechanisms behind the highaltitude adaptations of vertebrates have been widely explored in humans (Huerta-Sanchez et al., 2014), antelopes (Ge et al., 2013), pigs , dogs (Gou et al., 2014), yaks (Qiu et al., 2012), snub-nosed monkeys (Yu et al., 2016), horses (Hendrickson, 2013), ground tits Qu et al., 2013) and fish .
The Tibetan chicken (TC; Gallus gallus (Linnaeus, 1758)) is a unique breed native to the Qinghai-Tibet Plateau that is well adapted to hypoxic conditions. This is reflected in their having a higher hatchability than lowland breeds when eggs are incubated at high altitudes (Visschedijk, 1985). Recent studies have explored the molecular mechanisms underlying various aspects of the adaptations of TCs to hypoxia. These have revealed physiological differences in TC embryos lower lactate dehydrogenase activity in the heart (Jia et al., 2016), lower water loss (Zhang et al., 2008) and lower complex I activity in the brain (Bao et al., 2007) than lowland chicken embryos incubated under hypoxic conditions. Several studies have also linked single-gene variations to the high-altitude adaptations of TCs. For example,  reported that the allele frequency of a non-synonymous SNP of the EPAS1 gene correlated with altitude, while a polymorphism (m.9441G > A) in the ATP-6 gene has been found that may affect the specific functions of the ATP-6 enzyme and may relate to the high-altitude adaptations of TCs (Zhao et al., 2016). Liu et al. (2017) also found that the haplotype diversity of the COX-III gene in TC populations was markedly higher in TCs than in lowland chicken populations and that the mitochondrial membrane was more hydrophobic. A further study by Su et al. (2016) found that the haplotypes H4 and HE4 correlated negatively with high-altitude adaptations in the Tas2rs gene. Additionally, studies suggested that chickens with the A allele at m.10081A > G had a probability of the ability to adapt to hypoxia over 2.6 times larger than those with the G allele . Studies analyzing variation in the TC genome and transcriptome Zhang et al., 2016;Tang et al., 2017;Jiang et al., 2018) and in TC genome methylation Peng et al., 2018), have identified a series of pathways and genes related to hypoxia in TCs. However, despite these numerous studies having previously explored the adaptation mechanisms of TCs, the transcriptomic basis of these adaptations has not been well established.
Skeletal muscle fibers represent one of the most abundant cell types in mammals (Ohlendieck, 2011). They can modify their size and their metabolic or contractile properties in response to various stimuli, such as neural activity, mechanical stress, metabolic and hormonal influences and environmental factors (Chaillou, 2018). A decline in oxygen availability (hypoxia) has been proposed to induce metabolic adaptations and the loss of mass in skeletal muscle tissue (Chaillou, 2018). This has been supported by human studies, which have shown that Tibetans have unusually high muscle myoglobin contents (Gelfi et al., 2004).
In this study, high-throughput sequencing was used to generate transcriptomes from the pectoralis tissues of TCs and Beijing fatty chickens (BFCs) raised in Tibet, to investigate the differences in gene expression between these two breeds that show extreme differences in hypoxic adaptation. Our objective was to identify differentially expressed genes (DEGs) that might be involved in regulatory functions related to hypoxic adaptations in TCs.

Animals and management
All chickens were provided by the Institute of Animal Science and Veterinary Science of the Tibet Academy of Agricultural and Animal Husbandry Sciences of China. One hundred TC and lowland BFC eggs were incubated and reared in Tibet at an altitude of 3700 m.
At the start of the incubation period, the eggs were placed in an incubator with their large end on top, to ensure that the beaks of the chicks could easily enter the air chamber to breathe at hatching, which reduced the risk of death. For the first 18 days of incubation, the temperature was set at 37.8°C and the humidity at 60%. The eggs were turned by 90°once every two hours. After 18 days of incubation, the eggs were transferred to the hatching tray, and for days 19 to 21 the temperature was set at 37.5°C and the humidity at 70%.
After hatching, the chicks were placed on the floor in separate cages in a wellventilated, windowless room at a stocking density of 12 chickens/m 2 . The room was equipped with shaded lamps (2 m high) placed at 4-m intervals, such that each of the lamps illuminated 20 m 2 . The lamps were fitted with 40-W (20 lux) bulbs for the first week post-hatching, and 25-W (10 lux) bulbs from the second week post-hatching, with the bulbs being frequently checked for damage to prevent irregularities in light intensity. Illumination was maintained at 10 lux for 23 hours/day from hatching until three days of age to ensure feeding and drinking, after which it was reduced to 10 lux for 18 hours/day until seven days of age. It was subsequently shortened by 10 minutes per day until an illumination period of 8 hours/day was reached, and this was maintained until 20 weeks post-hatching. The illumination period was then increased by 0.5 hours/week until it reached 15 hours/day and remained unchanged thereafter.
Feed and water troughs were provided at 6-cm intervals and feed and water were provided ad libitum. All chickens were kept in the same rearing environment until 300 days of age, at which time muscle tissue samples were collected from four TCs and four BFCs (sampled individuals were randomly selected).
All experiments were conducted according to the guidelines for the care and use of experimental animals established by the Ministry of Science and Technology of the People's Republic of China. Chicken rearing and tissue sample collection procedures were approved by the Animal Welfare Committee of Sichuan Agricultural University (approval number: DKY-S20160906). All efforts were made to minimize suffering.

Chicken tissue sample collection and mRNA sequencing
Pectoralis muscle tissue samples for mRNA sequencing (RNA-seq) were collected from eight healthy hens of similar body weight at 300 days of age. After collection the samples were immediately placed in liquid nitrogen and stored at −80°C until later use. The extraction of total RNA using the TRIzol reagent (Invitrogen Life Technologies, Inc., Beijing, China) was performed according to the manufacturer's instructions. The quality and integrity of the RNA was determined using agarose gel electrophoresis. Total RNA was purified using an RNeasy mini kit (QIAGEN, Valencia, CA), and the purity of the RNA was monitored using a Nanodrop-2000c Nucleic Acid Protein Analyzer (Thermo Scientific, Shanghai, China). FastQC (http://www.bioinformatics.bbsrc.ac.uk/projects/fastqc) was used to get a visual overview of the read quality. RNA sequencing was performed using an Illumina (San Diego, CA, USA) HiSeq 2000 sequencing system. Reads were mapped against the reference genome Galgal5.0 from Ensembl (http://www.ensembl.org/index. html), using TopHat version 2.0.7 (http://ccb.jhu.edu/software/tophat/index.shtml). Cufflinks version 2.2 (Cufflinks, http://cole-trapnell-lab.github.io/cufflinks/install/) was applied to quantify gene expression and obtain transcripts per million (TPM) expression values (Trapnell et al., 2010). Abundance files were generated by applying Cuffquant (Cufflinks, RRID: SCR014597) to the read mapping results. Log 2-transformations of the TPM + 1 values were used for subsequent analyses for genes with >10 TPM in more than 80% of the samples (Mortazavi et al., 2008). The raw reads were submitted to the National Center for Biotechnology Information (NCBI, https://www.ncbi.nlm.nih.gov/).

Statistical analysis
We evaluated all DEGs using the following steps. DEGs were identified and cluster analysis of the two groups was performed using MeV (version 4.9.0, http:// mev.tm4.org/) (Saeed et al., 2006). We used t-tests to obtain DEGs, with adjusted P -values 0.01 and | log 2 (Fold Change)| 1 being considered statistically significant. Given the small sample size, we used a relatively lenient cut-off to guard against false negative errors. The online biological information annotation database (http://metascape.org) was used for functional annotation of the DEGs (Zhou et al., 2019). Gene Ontology (GO, http://geneontology.org/) and Kyoto Encyclopedia of Genes and Genomes (KEGG, https://www.genome.jp/kegg/) pathway enrichment analyses were performed for all DEGs (Kanehisa et al., 2008).

RNA quality detection and reference sequence alignment
Sequencing experiments need high-quality samples, with an RNA integrity number (RIN) greater than 7 being required for Illumina sequencing. The RINs of the samples used in this study were all >7, and the RNA of the eight sample tissues met all the requirements for sequencing. The quality analysis of the transcriptome sequencing data from the chicken pectoralis muscle tissues is shown in table 1. We generated 80.72 Gb of RNA-seq data from eight samples: four TCs and four BFCs, all of which were raised at a high altitude in Tibet. More than 80% of the reads in each sample could be reliably aligned to the chicken reference genome (data not shown). The expression data from the current study have been submitted to the Gene Expression Omnibus (GEO, http://www.ncbi.nlm.nih.gov/geo).

Gene expression and differentially expressed gene (DEG) analysis
TPM values were used to estimate the expression levels of the identified genes (Mortazavi et al., 2008). We found that, on average, 67.22% (12,333) of the annotated protein-coding genes in the eight samples had TPM expression values greater than 10. For the analysis of up-or downregulated genes, based on a conservative | log 2 (FC)| 1 threshold and similar biological alteration patterns in the two chicken breeds under hypoxic conditions, 48 DEGs were obtained ( fig. 1, table 2). Genes with higher expressions in TCs than in BFCs were defined as upregulated genes; genes with lower expressions in TCs were defined as down-regulated genes. Under the hypoxic incubation conditions used in this study, all DEGs were upregulated and had higher expression levels in TCs (table 2).

Figure 1.
Volcano plot of genes from four Tibetan chickens and four Beijing fatty chickens, reared at a high altitude in Tibet. Red dots represent genes that were differentially expressed genes, being significantly higher in TC; blue dots indicate that the genetic differences were not significant (and no down-regulation was observed).

Gene ontology and KEGG pathway enrichment analysis
The main gene enrichment analyses presented here ( fig. 2) were conducted on the Metascape website. Consequently, each group of GO terms has a similar biological meaning, as they share similar gene members. The higher-ranked annotation groups most likely have consistently lower P -values for their annotation members. Interestingly, these genes are mainly associated with immune system responses and ion channel activity in response to stimuli. P -values were used for gene enrichment analysis, with a smaller P -value indicating more significant biological enrichment. In the classification of enrichment, the biological processes enriched included: branching morphogenesis of epithelial tubes, heart morphogenesis, and skeletal system development. These enriched entries indicate that to maintain their metabolism, growth and development, TC organization continues to synthesize and metabolize protein complexes, support cell membrane synthesis and skeletal muscle development, and generate intermediate products from various energy metabolism processes. They also resist and respond to various harmful substances produced by the hypoxic stimulation. KEGG pathway enrichment analysis showed no significantly enriched pathways associated with the identified DEGs. Hypoxic adaptation requires the interaction of many biochemical pathways and synergistic changes in gene expression and structure. Therefore, DEGs in KEGG pathways may have important physiological functions in biological processes. Only nine DEGs were annotated in KEGG pathways, as listed in table 3. Among these, GRK5 and DNM3 are involved in the endocytosis KEGG pathway, CACNA1H participates in the MAPK and calcium signaling pathways, and TGIF2 plays a role in the TGF-beta signaling pathway (table 3). These pathways are closely related to the hypoxic adaptation process, and these genes may therefore be candidate hypoxia-related genes with important physiological functions.

Discussion
TCs are uniquely adapted to the extreme high-altitude environment that they inhabit . A previous comparative genomics study showed that genes positively selected for in highland chicken populations were related to cardiovascular and respiratory system development, which indicates strong adaptation to conditions of hypoxia and high-intensity solar radiation . These authors computed mean expression values between highland and lowland chickens to explore the unique adaptations of TCs to hypoxia and high-dose ultraviolet radiation in high-altitude environments, using two biological replicates . Here, we explored TCs' adaptations from a muscle transcriptome perspective, and used four TCs and four BFCs in each group -a sample size that was sufficient for analysis. We identified gene expression changes related to TCs' hypoxic adaptation mechanisms and comprehensively investigated the pathways in which DEGs were involved. During hypoxia, an organism maintains blood oxygen levels through systematic regulation mechanisms. At the molecular level, cells in a hypoxic environment may sense hypoxic signals through oxygen receptor molecules on their membranes, and in response, induce hypoxic response genes that may further alter the expression of other genes, including cytokines, lymphokines, and hormones.
Balancing the oxygen supply and demand is necessary for cell metabolism and survival, and animals use various mechanisms to maintain this balance in the body under hypoxic conditions. For example, erythropoietin and vascular endothelial growth factor are induced to activate increased oxygen delivery, thus compensating for the reduced oxygen supply (Kimura et al., 2000;Figueroa et al., 2002). Another major intracellular hypoxia adaptation is to switch from using oxidative phosphorylation to glycolysis as the primary way to obtain ATP. Cells can thus respond to hypoxia both internally and externally, by altering their gene expression to maintain homeostasis. Indeed, it has been further reported that under hypoxic conditions, the expression of genes encoding enzymatic components of the respiratory chain is decreased, while the expression of genes encoding glycolytic components is increased (Semenza et al., 1996).
We have previously shown that the EPAS1, ATP-6, COX-3, Tas2rs and mt-CO3 genes are related to hypoxia, and that there are SNPS in these genes associated with plateau adaptation Su et al., 2016;Zhao et al., 2015); however, these genes were not differentially expressed in our current study. Instead, our results showed differential expression of B4GALT5, GRK5, wisp2, EIF2AK1 and CDC45 in the breast muscles of TCs and BFCs, and we thus speculate that these genes might also be related to hypoxic conditions as their expression levels in the pectoralis muscles of TCs were significantly upregulated under hypoxic conditions. The WNT-1 inducible signaling pathway protein 2 (WISP-2), encoded by the wisp2 gene, is a secreted protein member of the connective tissue growth factor/cysteine-rich 61/nephroblastoma overexpressed (CCN) family (Fuady et al., 2014). Wisp2 is regulated exclusively by hypoxia-inducible factor-2 (HIF-2) (Stiehl et al., 2012), and has been found to be expressed in adult skeletal muscle, colon, ovary tissues, and in the fetal lung (Pennica et al., 1998). Its expression correlates inversely with the aggressiveness of breast, pancreatic, and colon cancer, which suggests it has tumor suppressor-like activity (Banerjee and Banerjee, 2012). The upregulation of wisp2 that we observed here may be related to the TCs' adaptation to hypoxia. The expression level of the EIF2AK1 (eukaryotic translation initiation factor 2-alpha kinase 1) gene in U87 glioma cells is downregulated in response to hypoxic conditions (Minchenko et al., 2016). In this study, EIF2AK1 was more highly expressed in TCs (TPM: 173.50) than in BFCs (TPM: 81.00). A previous study suggested that GRK5 (G protein-coupled receptor kinase 5) deficiency makes mice more susceptible to various behavioral and cognitive impairments (Singh et al., 2016). Further research is needed to confirm whether the EIF2AK1 and GRK5 genes are associated with the high-altitude adaptations of TC.
B4GALT5, GRK5, wisp2, EIF2AK1, and CDC45 expression levels in the pectoralis muscles of TCs were significantly increased compared with those in BFCs under conditions of high altitude and low oxygen. These results suggest that hypoxia can upregulate the expression of several kinase and phosphatase genes in the muscle tissues of TCs (Minchenko et al., 2016;Singh et al., 2016), thereby inducing the synthesis of proteases, promoting cell development, and producing sufficient amounts of ATP through protein synthesis and metabolism to meet the needs of energy metabolism. This cellular response to environmental stimuli is one of the ways in which TCs' adaptability to high altitude hypoxic environments can be regulated at the molecular level.
Major histocompatibility complex (MHC) class I molecules are cell-surface receptors that present peptides from intracellular antigens to cytotoxic T lymphocytes, and serve as recognition elements for natural killer (NK) cells. MHC class I polygeny and polymorphism contribute to the breadth of the immune response to pathogens. Chickens possess two MHC gene clusters, the B locus and the Rfp-Y locus. The chicken B locus, called the "minimal essential MHC" because it carries all the hallmarks of the mammalian MHC, is largely contained within a 92-kb region (Kaufman et al., 1999). Two MHC class I genes flank either side of TAP1 and TAP2, while genes encoding LMP2 and LMP7 are absent from the MHC region, and may have been lost from the genome entirely (Kaufman et al., 1999). The proximity of TAP to class I genes may allow for the co-evolution of these proteins (Kaufman et al., 2015), and it has been suggested that TAP transports peptides of the appropriate specificity for presentation on the linked MHC class I molecules. Consequently, the upregulation of TAP2 in muscle tissue from TCs (TPM: 21.5) relative to BFCs (TPM: 8.25) found in this study may facilitate respiration during hypoxia.
During the development of avian muscle, there is serial expression of myogenin (MYOG) isoforms within the muscle fiber. In this study, the expression level of MYOG was elevated in TCs (TPM: 76.5, relative to 36 for BFCs) under hypoxia. This increased level of MYOG may suggest that alterations in MYOG isoform transitions were induced by hypoxia. Furthermore, it has been reported that activation of satellite cells may be responsible for the upregulated expression of MYOG, which stimulates compensatory muscular hypertrophy in animals (Snow, 1990). This could explain why lowland chicken breeds tend to exhibit hypertrophy when reared at high altitudes, whereas TCs continue to develop normally when exposed to chronic hypoxia.
While screening the DEGs in TCs and BFCs, some of the genes were associated with GO entries involved in nephron morphogenesis and nephron tubule development in response to external stimuli. Of the 42 enriched GO entries, half were related to renal tubule development and morphogenesis, which highlights the possible effects of hypoxia on kidney function. This effect of tubular hypoxia on the proximal tubule may be linked to the high energy requirements and reliance on aerobic metabolism of this structure (Gilbert, 2017). A previous study demonstrated that when SGLT was inhibited, renal cortex Po2 in the diabetic rat kidney was normalized, which implies that increased proximal tubule transport contributes to the development of hypoxia in the diabetic kidney (O'Neil et al., 2015). Chronic hypoxia is also an important factor contributing to the extent of glomerular injury in the tubulointerstitium (Theilig, 2010).
In summary, we obtained the transcriptome profiles of pectoral muscle tissues from TCs and BFCs under hypoxic conditions in Tibet. Our results indicate that hypoxia is related to renal morphogenesis and renal tubule development. Only nine of the DEGs identified were linked to KEGG pathways, including the endocytosis pathway, and the MAPK, calcium, and TGF-beta signaling pathways. Notably, Zhang et al. (2016), using genomic selection signature and transcriptome profile analyses to identify candidate genes involved in TCs' adaptation to highland environments, also found that both positive-selection genes and DEGs were clustered in the MAPK signaling pathway . Although much work still needs to be done, our study provides a description and interpretation of the functions and regulation of genes relevant to TCs' hypoxic adaptations.