Comparative Transcriptome Analyses of Gayal (Bos frontalis), Yak (Bos grunniens), and Cattle (Bos taurus) Reveal the High-Altitude Adaptation

Gayal and yak are well adapted to their local high-altitude environments, yet the transcriptional regulation difference of the plateau environment among them remains obscure. Herein, cross-tissue and cross-species comparative transcriptome analyses were performed for the six hypoxia-sensitive tissues from gayal, yak, and cattle. Gene expression profiles for all single-copy orthologous genes showed tissue-specific expression patterns. By differential expression analysis, we identified 3,020 and 1,995 differentially expressed genes (DEGs) in at least one tissue of gayal vs. cattle and yak vs. cattle, respectively. Notably, we found that the adaptability of the gayal to the alpine canyon environment is highly similar to the yak living in the Qinghai-Tibet Plateau, such as promoting red blood cell development, angiogenesis, reducing blood coagulation, immune system activation, and energy metabolism shifts from fatty acid β-oxidation to glycolysis. By further analyzing the common and unique DEGs in the six tissues, we also found that numerous expressed regulatory genes related to these functions are unique in the gayal and yak, which may play important roles in adapting to the corresponding high-altitude environment. Combined with WGCNA analysis, we found that UQCRC1 and COX5A are the shared differentially expressed hub genes related to the energy supply of myocardial contraction in the heart-related modules of gayal and yak, and CAPS is a shared differential hub gene among the hub genes of the lung-related module, which is related to pulmonary artery smooth muscle contraction. Additionally, EDN3 is the unique differentially expressed hub gene related to the tracheal epithelium and pulmonary vasoconstriction in the lung of gayal. CHRM2 is a unique differentially expressed hub gene that was identified in the heart of yak, which has an important role in the autonomous regulation of the heart. These results provide a basis for further understanding the complex transcriptome expression pattern and the regulatory mechanism of high-altitude domestication of gayal and yak.


INTRODUCTION
Species living at high altitudes are exposed to strict selection pressures and physiological challenges owing to harsh environmental conditions, such as thin air, cold temperatures, ultraviolet exposure, and low pressure (Miao et al., 2015).Despite the harsh conditions, numerous species including Tibetan pigs (Jia et al., 2016;Zhang B. et al., 2017), Tibetan sheep (Zhang et al., 2013), Tibetan chickens (Gou et al., 2014), and yak (Qiu et al., 2012;Qi et al., 2019) have evolved unique physiological characteristics, such as superior blood oxygen transport system and high metabolic efficiency, to adapt to the harsh living pressure (Lan et al., 2018), as have the native Tibetan people (Ge et al., 2012;Tashi et al., 2014).At present, published studies have identified EPAS1, EGLN1, and PPARA genes that play important roles in high-altitude adaptation (Haasl and Payseur, 2016;Heinrich et al., 2019).Exploring the molecular mechanisms underlying hypoxia adaptation has long garnered attention.
Gayal, also known as Drung cattle (Bos frontalis) (Winter et al., 1986;Uzzaman et al., 2014), is a unique semi-wild cattle breed, mainly found in the typical alpine valleys and subtropical rain forests in the Drung and Nujiang river basins of Yunnan Province, southwestern China.Studies have shown that the level of hemoglobin in gayal, as well as the number of red blood cells and white blood cells, are equivalent to those of yak, which is the physiological basis for the ability of gayal to resist invasion by bacteria, viruses, and parasites in the alpine environment (Tian et al., 1998).Yak (Bos grunniens), a native mammal on the Qinghai-Tibet Plateau and its adjacent regions, provides meat and other necessities for Tibetans.Compared with lowland cattle, yaks have a larger alveolar area per unit area, thinner alveolar spacing, thinner gas-blood barrier (Wei and Yu, 2008), and stronger expression of some genes related to oxygen supply as well as defense system under hypoxia pressure (Wang et al., 2016).Therefore, a thorough understanding of the genetic basis of the physiological characteristics of gayal and yak will provide insight into their adaptation to the high-altitude environment.
With recent rapid progress in next-generation genome sequencing (NGS) technologies, high-throughput RNA-Sequencing (RNA-seq) technology represents a powerful and cost-efficient approach to explore species domestication, breeding, and genetic molecular mechanisms in organisms (Salleh et al., 2017;Keren et al., 2018).RNA-seq can detect and quantify gene expression with digital measurements, which are especially sensitive for low-expressed genes (Mortazavi et al., 2008).In addition, it is also used to improve gene annotation, discover novel genes or transcripts, and survey alternative splicing (AS) events at single-nucleotide resolution (Sultan et al., 2008;Ozsolak and Milos, 2011).Thus far, transcriptome studies on gayal and yak have been widely conducted.The muscle transcriptome analysis of Indian Mithun showed that hub genes including MTMR3, CUX1, LONRF3, PLXNB2, KMT2C, ZRSR2Y, PRR14, USP9Y, SLMAP, and KANSL2 might contribute to muscle growth (Mukherjee et al., 2020).Transcriptomic analysis of yaks living at different altitudes indicates PI3K-Akt, HIF-1, focal adhesion, and ECM-receptor interaction pathway were significantly enriched, and the EPAS1 expression increases with altitude (Qi et al., 2019).Comparative transcriptome analysis of four organs (heart, kidney, liver, and lung) between yak and cattle revealed that DEGs associated with the blood supply system, modulation of cardiac contractility, vascular smooth muscle proliferation, and the glutamate receptor system probably play an important role in yak adaptation to hypoxia environments (Wang et al., 2016).However, gayal is a species that is well adapted to their local alpine valley environment, and few transcriptomic studies focusing on regulation analysis of gayal in adaptation to the local high-altitude environment have so far not been reported.In addition, studies have found specific adaptive genes and mechanisms are distinct between species and populations (Heinrich et al., 2019;Pamenter et al., 2020).Therefore, we seek to understand how changes in altitude affect the transcriptional regulation of gayal, and whether the relevant transcriptional regulation mechanism is the same as that of yak, which shares a bovine subfamily.
In this study, we characterized gene expression profiles using RNA-seq data from six major tissues including heart, lung, liver, kidney, spleen, and muscle among gayal, yak, and cattle.We performed a comparative transcriptome analysis of six tissues in gayal vs. cattle and yak vs. cattle, respectively, and identified candidate DEGs related to adaptability to the plateau environment, as well as analyzed their functions and expression patterns.In addition, weighted gene co-expression network analysis (WGCNA) was further used to uncover the hub genes that regulate tissue function in high-altitude adaptation.The results could lay a foundation for further explorations of the genetic changes underlying hilly adaptations in subtropical mammals.

Sample Collection
Six tissues (heart, kidney, liver, lung, skeletal muscle, and spleen) were collected from three adult female gayal living in the semiwild preserved field at an altitude of 2,300-3,500 m in the Drung and Nujiang river basins.All tissue samples were taken and immediately snap-frozen in liquid nitrogen.In addition, the corresponding tissue transcriptome data of yak and cattle were retrieved from a previous publication (Tang et al., 2017) from the GEO databases under the accession numbers GSE93855 and GSE77020 (Supplementary Table S1).In addition, one gayal multi-tissue transcriptome data that we previously published was also included in this study to expand the sample size, which was stored in National Genomics Data Center under the accession code PRJCA002143.An overview of the samples used and the associated statistics are provided in Supplementary Table S1.

RNA Extraction, Library Construction, and Transcriptome Sequencing
Total RNA was isolated and purified with the TRIzol reagent (Life Technologies, Carlsbad, United States) according to the manufacturer's protocols.RNA purity was qualified using NanoDrop ND 2000 spectrophotometer at 260 and 280 nm (Thermo Fisher Scientific, Wilmington, United States), and RNA integrity was assessed using Agilent 2,100 Bioanalyzer (Agilent Technologies, Palo Alto, United States).The OD260/ 280 ratios of all samples were greater than 1.8, and the RNA integrity number (RIN) value of >7.0.Then, mRNA was broken down into small fragments using a magnesium RNA cleavage module (NEB, Ipswich, United States).Random primers and related reverse transcriptase were used to synthesize the firststrand complementary DNA (cDNA).Subsequently, the secondstrand cDNA was synthesized.The average insert size for the final cDNA library was 300 ± 50 bp.Finally, we conducted the 2 × 150 bp paired-end sequencing (PE150) on an Illumina Novaseq ™ 6,000 (LC-Bio Technology CO., Ltd., Hangzhou, China) following the supplier's recommended protocol.

Identification of Orthologous Genes
To identify orthologous genes among gayal, yak and cattle, we first downloaded protein sequence data of yak and cattle from the Ensemble genome browser 104 (http://asia.ensembl.org/index.html) and combined with the protein sequence data of the gayal reference genome assembled by ourselves (the accession code of National Genomics Data Center is PRJCA004132) for downstream analysis.We used ORTHOFINDER v2.2.6 (Emms and Kelly, 2015) software to delimit orthologous groups, in which all predicted protein sequences were compared using a BLAST all-against-all search (Camacho et al., 2009).The single-copy genes, duplicated genes, and species-specific genes were extracted from the ORTHOFINDER output, respectively.Among them, the 1:1:1 homologous genes of the three species were defined as the conserved single-copy homologous genes of each species, which were retained for the following analysis.

Data Quality Control, Processing, and Normalization of Gene Expression Levels
To ensure the accuracy of subsequent biological information analysis, the quality of the raw data was first checked using FastQC (http://www.bioinformatics.babraham.ac.uk/projects/ fastqc/).According to the evaluation results of raw reads, the software FASTP (https://github.com/OpenGene/fastp)with default parameters was used to remove adapters, reads with a high proportion of unknown nucleotides, and low-quality sequences (≤Q20) to obtain high-quality clean reads.To avoid biases due to sequence divergence across bovine subfamily genomes, we downloaded yak (BosGru_v2.0,http://ftp.ensembl.org/pub/release-104/fasta/bos_mutus/)and cattle (ARS_UCD1_2, http://ftp.ensem-bl.org/pub/release-104/fasta/ bos_taurus/) reference genomes from Ensemble database, and retrieved gayal (Drung_v1.2) reference genomes from National Genomics Data Center (NGDC) (https://bigd.big.ac.cn/).Then, the clean reads for each sample were mapped to the corresponding reference genome using HISAT2 software (Pertea et al., 2016).FeatureCount from the Subread package (Version 2.0.0)(Liao et al., 2014) was applied to calculate the numbers of reads mapped to each gene.Read pairs were counted only if they were uniquely mapped and properly paired.The gene expression level estimates may be biased among species due to factors such as the size of mRNA transcript and the qualities of different genome annotations (Blake et al., 2020).To avoid these issues, we used a custom script to only retain reads that mapped to the 1:1 orthologous which can be used for each of the three genomes (Blake et al., 2020).Taking into account the difference in sequencing depth of samples from different sequencing platforms, we first filtered out genes with low expression, and only keep genes with expression level CPM> 1 in at least half of the samples (Blake et al., 2020).We normalized the raw read counts using the weighted trimmed mean of M-values algorithm (TMM) in the edgeR package (Robinson et al., 2010).Then scaling factors were used to scale the expression levels of all orthologous.The normalized expression data for the trimmed orthologous were used for subsequent gene expression analysis.

Gene Expression Profile Analysis
After performing normalization, we obtained the TMM-normalized log 2 -transformed CPM values for each gene.The principal component analysis (PCA) was performed with the prcomp function in R using normalized log 2 -transformed CPM values.The hierarchical clustering of Pearson's correlation coefficients for each pair of samples was performed using the complete-linkage agglomerative method on the correlation distance matrix and generated symmetrical heatmap plot using the R gplots package.

Differential Expression Pattern Analysis
To characterize the high-altitude adaptive gene expression pattern of gayal and yak, we performed differential expressed analyses among heart, kidney, liver, lung, skeletal muscle, and spleen tissues in gayal vs. cattle group and yak vs. cattle group by using R package edgeR (Robinson et al., 2010).Genes with Benjamini Hochberg false discovery rate (FDR) values less than 0.05 and fold change more than 2 were identified as DEGs.If the log 2 fold changes more than 1 were defined as up-regulated DEGs, or less than -1 were defined as down-regulated.We performed Gene Ontology (GO) function annotation and KEGG pathway enrichment analysis for DEGs by the Database for Annotation, Visualization and Integrated Discovery (DAVID) v6.8 using the Bos taurus genome as the background set (Huang Da et al., 2009).The significantly enriched GO and KEGG terms using a 0.05 cutoff for the p-value with Benjamini-Hochberg test adjustment.The jvenn and UpSetR R package (version 1.4.0)were used to identify and visualize the DEGs shared by two comparison groups across the multiple tissues and between two comparison groups in the same tissue, respectively (Conway et al., 2017).In addition, to better understand the hypoxia adaptation of species, we also collected a set of known or potential high-altitude adaptation candidate genes from previously published literature (Zhang et al., 2014;Wang et al., 2016;Peng et al., 2017).

Construction of Gene Co-Expression Network
The adaptation of species to the hypoxia environment is a complex biological process, which involves the coordinated regulation of multiple genes.WGCNA is a feasible method to clarify the regulatory relationship between genes and identify important hub genes based on gene expression profiles to explore the mechanism of hypoxia adaptation (Wang et al., 2021).Therefore, to further identify the relevant regulatory gene modules related to the functions of each tissue and explore the core driver genes, we used the CPM values of filtered genes to perform WGCNA analysis on three species, respectively (Langfelder and Horvath, 2008).The hierarchical cluster analysis on all tissues was carried out using the hclust function.The soft-power threshold beta was determined by the function "sft$powerEstimate" (Chen et al., 2019).The minimal gene module size was set to 30 to obtain appropriate modules, and the threshold to merge similar modules was set to 0.25.Gene modules were detected based on the TOM matrix.The correlations between modules eigengene and tissues were investigated using the Pearson correlation, and only modules with correlation coefficients >0.7 and p < 0.001 were considered tissue-relevant modules (Touma et al., 2016).Here, we designate the highest top 5% of genes as the hub genes based on the intramodular connectivity in each functional module, which is calculated by summing the connection strengths with other intramodular genes (Yang et al., 2014).Functional enrichment analyses for tissue-related module genes were performed using DAVID v6.8 using the Bos taurus genome as the background set, and the p-values were corrected for multiple testing using Benjamini-Hochberg correction (Huang Da et al., 2009).

Summary Statistics of RNA-Seq Data
The transcriptome data of six tissues (heart, kidney, liver, lung, skeletal muscle, and spleen) from four gayals, three yaks, and three cattle were used in this study (Supplementary Table S1).One gayal heart sample (gayal3_heart) was excluded due to low sequencing quality.In total, ~327.55G clean data has been generated after removing adapter-embedded or low-quality reads.We calculated Q20 (99% base call accuracy), Q30 (99.9% base call accuracy), and GC-content of the clean data.The average values of Q20 and Q30 in the samples were 97.04% (ranging from 94.60 to 98.47%) and 92.38% (ranging from 89.39 to 95.43%), respectively (Supplementary Table S2).The average base GC content was 50.80%, ranging from 47.69 to 55.03%, indicating that the base composition and quality of the sequencing data were qualified.Then the high-quality reads of each sample were mapped to the corresponding reference genome (i.e., Drung_v1.2,BosGru_v2.0,ARS_UCD1_2) using the HISAT2 software (Pertea et al., 2016).The average reads mapping rates in gayal, yak, and cattle were 86.35, 93.60, and 95.58%, respectively.Detailed statistics and mapping information for each sample were summarized in Supplementary Table S2.

Ortholog Identification and Gene Expression Levels Normalization
To reduce the effect of chromosome number difference (Bos frontalis, 2n = 58; and Bos grunniens, 2n = 60; Bos taurus, 2n = 60) and the qualities of different genome annotations among three bovine species, we first used the blast function of ORTHOFINDER2 (Emms and Kelly, 2015) to identify orthologous genes among the gayal, yak and cattle.A total of 19,174 protein orthogroups were identified.Among them, only 6,371 orthogroups were high-confidence single-copy orthogroups of three species (Supplementary Table S3).Then the corresponding gene names were extracted by a custom script and retained reads that mapped to the single-copy homologous genes of three species genome, and the ensemble gene names of cattle were uniformly used for subsequent analysis.In the 6,371 single-copy orthologous gene expression datasets, genes with low expression (CPM <1) in more than half of the samples were further excluded, and 4,712 genes were obtained (Blake et al., 2020).Then we normalized 4,712 single-copy orthologous genes using the weighted TMM of edgeR (Hao et al., 2019).

Gene Expression Profiles Across Tissues and Bovine Species
To investigate broad gene expression patterns in multiple tissues across different bovine species, we performed hierarchical cluster analysis based on the gene expression of six tissues among gayal, yak, and cattle.As shown in Figure 1A, the tissue-specific expression pattern of lung, liver, spleen, and kidney shows that the same tissues of different bovine species gather earlier than different tissues of the same bovine species.The result indicated that tissue differentiation may precede species differentiation, which was consistent with previous studies (Necsulea and Kaessmann, 2014;Hao et al., 2019).However, the samples of heart and muscle from each species were more similar than the same tissue between different species (speciesspecific expression pattern), which reflects striated muscular tissues may have similar expression profiles (Merkin et al., 2012;Tang et al., 2017;Feng et al., 2020) (Figure 1A).In addition, our expression profiles suggested that the differences in global gene expression among tissues were more significant than that among altitudes or the same tissue of different species, which was consistent with the results of previous studies (Necsulea and Kaessmann, 2014;Tang et al., 2017;Hao et al., 2019).Furthermore, principal component analysis (PCA) was performed on the gene expression of 59 samples.The score plot of PCA analysis for six main organs shows that 45.4% of the variance can be explained by the first two principal components, accounting for 26.4 and 19% of the variance, respectively (Figure 1B).In essence, 59 samples were clustered based on similar tissue types, which confirmed that tissue differentiation may precede species differentiation.

Gene Expression Changes Accompanying High-Altitude Environment
To detect expression changes of gayal and yak in response to the high-altitude stress, we used edgeR (Robinson et al., 2010) to identify DEGs within each tissue.According to the tissue and species we considered, 1,218 to 1,585 interspecies DEGs were identified in the comparison of six type tissues in gayal and cattle (Table 1).Among the six tissues, the heart has the largest number of DEGs (618 up-regulated and 967 down-regulated genes).By upset analysis of all DEGs in each tissue, 446 DEGs shared in six tissues were identified, and each gene showed consistent upregulation or down-regulation in all tissues (Figure 2A).In the yak vs. cattle group, we identified 1,995 DEGs among six types of tissues.The number of DEGs among tissues ranges from 543 in liver tissue (202 up-regulated and 341 down-regulated genes) to 919 in lung tissue (387 up-regulated and 532 down-regulated genes) (Table 1).There were 154 DEGs with consistent regulation direction shared among the six tissues, of which 64 genes were also present in all tissues of the gayal vs. cattle group (Figure 2B).
To further investigate the potential biological function of DEGs for environmental adaptation in gayal and yak, we performed functional annotation analysis based on the DAVID database.In the gayal vs. cattle group, DEGs were annotated in energy metabolism, cardiovascular system, blood coagulation system, and immune system among tissues (Supplementary Table S4).As such, the DEGs of the heart were enriched in the oxidation-reduction process (GO: 0055114), positive regulation of angiogenesis (GO:0045766), positive regulation of blood coagulation (GO:0030194).Muscle DEGs were enriched in blood coagulation (GO:0007596), patterning of blood vessels (GO:0001569), skeletal muscle cell differentiation (GO:0035914) (Figure 3).We also observed DEGs in the gayal vs. cattle group were significantly enriched in hypoxia response, metabolism, and immune system pathway.For instance, hematopoietic cell lineage (bta04640) was enriched in four tissues (heart, spleen, lung, and kidney) and galactose metabolism (bta00052) was enriched in four tissues (heart, liver, spleen, and lung) (Supplementary Table S4).In the yak vs. cattle group, DEGs identified across six tissues pairwise comparisons were also enriched in the cardiovascular system, energy metabolism, and blood coagulation system, which included response to oxidative stress (GO:0006979) and oxidation-reduction process (GO:0055114) in the heart, and patterning of blood vessels (GO:0001569) and positive regulation of canonical wnt signaling pathway (GO:0090263) in the muscle (Figure 4).Similar results were found in the KEGG analysis, pathways related to energy metabolism and the immune system were significantly enriched (Supplementary Table S4).
To explore the similarities and differences of expression regulation genes related to environmental adaptability among gayal and yak, we further analyzed the shared and unique DEGs  between the gayal vs. cattle group and the yak vs. cattle group.Among the six tissues of gayal and yak, the heart and lung shared the largest number of DEGs, with 397 (101 up-regulated and 296 down-regulated genes) and 396 genes (146 up-regulated and 250 down-regulated genes), respectively (Figure 5).In addition, we have identified 822-1,331 and 298-500 unique DEGs in the  the intrinsic apoptotic signaling pathway in response to DNA damage (GO:0008630).Liver-shared DEGs were enriched in innate immune response (GO:0045087).Spleen-shared DEGs were enriched in chemotaxis (GO:0006935).And the shared DEGs in the muscle were enriched in skeletal muscle cell differentiation (GO:0035914) (Supplementary Table S4).As for the unique DEGs of the gayal vs. cattle group, we found that five unique DEGs in the heart related to apoptosis (bta04210), include CASP9, ENDOG, AIFM1, CASP6, TNFSF10, RIPK1, FASLG, and TNF.Seven unique genes (BLOC1S6, SERPINE2, TFPI2, HPS4, SERPING1, GP1BA, and TFPI) are associated with blood coagulation (GO:0007596) in the lung.Spleen-unique DEGs were enriched in lymphocyte chemotaxis (GO:0048247) and inflammatory response (GO: 0006954).Furthermore, we found that positive regulation of angiogenesis (GO:0045766) was significantly enriched in specific genes in the lung, kidney, and spleen; hematopoietic cell lineage (bta04640) was also significantly enriched in the unique genes of the heart, spleen, lung, and kidney.In addition, we found that the unique DEGs in the heart, liver, and spleen were also related to galactose metabolism (bta00052).These unique DEGs may play an important role in the adaptation of gayal to the alpine and valley environment (Supplementary Table S4).Among the unique DEGs of the yak vs. cattle group, we found that four genes (CD40LG, SRF, F2R, and F5) related to platelet activation (GO:0030168) in the heart.Blood coagulation (GO:0007596) related genes (PROCR, THBD, F10, PROS1, F3, and KNG1) were also significantly enriched in the lung.The unique genes in the spleen were enriched in immune-related items, such as the NOD-like receptor signaling pathway (bta04621), and TNF signaling pathway (bta04668).In addition, cell adhesion molecules (CAMs) (bta04514) were uniquely enriched in the heart and lung of yak, which was reported as the intensity of hypoxia increases lead to marked reduced cell adhesion (Kaiser et al., 2012) (Supplementary Table S4).

Expression Regulation of High-Altitude Related Genes in the DEGs
To further establish the hypoxia adaptation of gayal and yak, we collected a candidate gene set with known or potential functions in hypoxia responses from previously published literature (Zhang et al., 2014;Peng et al., 2017).Of those, 295 and 187 genes were DEGs in at least one tissue of gayal and yak, respectively (Supplementary Table S5).Among which the heart has the largest number of high-altitude related genes in both comparison groups DEGs (149 in gayal vs. cattle group and 85 genes in yak vs. cattle group).Notably, in the altitude-related gene set, we found that several positively selected genes related to hypoxia response in Tibetans were differentially expressed in two comparison groups.For example, EGLN1, also known as PHD2, as a known target gene of the hypoxia-sensing pathway (Haasl and Payseur, 2016;Heinrich et al., 2019;Storz, 2021), showed a consistent down-regulation across all tissues of yak and gayal compared with lowland cattle, but the expression level and pattern vary in different tissues.HFE (Homeostatic Iron Regulator), a gene involved in indirectly regulating iron balance in animals via reducing transferrin-mediated iron uptake (Salter-Cid et al., 1999), was uniquely down-regulated in all tissues of gayal except for kidney compared with the lowland cattle.Additionally, HOXB6, a gene associated with the regulation of the erythropoietic system, was particularly up-regulated in the spleen and kidney of the yak (Zimmermann and Rich, 1997).

Tissue-specific Module Detection and Hub Genes Identification
To explore genes and co-expression modules with similar expression profiles related to high-altitude adaptation in various tissues, WGCNA analysis was independently performed on the three species based on the CPM values of 4,712 filtered genes.The sample clustering dendrogram of the gayal, yak, and cattle are shown in Figures 6A, 7A, 8A, respectively.In both species, all samples were included in the clusters.The soft-power threshold of 10, 18, and 16 (scale-free R 2 of 0.85) were selected for further analysis of the gayal, yak, and cattle, respectively (Figures 6B, 7B, 8B).Next, the gene modules were detected based on the topological overlap matrix (TOM), and sixteen modules were detected under the gayal, while ten modules were detected under the yak (Figures 6C, 7C, 8C).To identify tissue-specific modules, we calculated the correlation coefficients between the module and each tissue (Figures 6D, 7D, 8D).Further, GO and KEGG analyses were conducted to investigate the biological function of the genes in each tissuerelated module.
Indeed, the relevant modules enriched in six tissues are highly related to the tissue-specific functions.For example, in the tissuerelated modules of gayal, the MEcyan module was significantly correlated with heart (r = 0.88, p = 3 × 10 -8 ), which was enriched in functional annotations such as negative regulation of cardiac muscle cell apoptotic process (GO:0010667), and pathway in MEcyan module were involved in cardiac muscle contraction (bta04260).Lung-related module (MEgreenyellow, r = 0.99, p = 4 × 10 -18 ) was enriched in HIF-1 signaling pathway (bta04066) (Supplementary Table S6).Similar results were found in the yak, MEgreenyellow was the most significantly correlated with the heart of yak (r = 0.98, p = 3 × 10 -12 ), which was enriched in functional annotations such as positive regulation of cardiac muscle cell proliferation (GO:0060045), positive regulation of osteoblast differentiation (GO:0045669), and KEGG including Adrenergic signaling in cardiomyocytes (bta04261).The module of MEbrown was significantly correlated with the lung of yak (r = 0.99, p = 3 × 10 -14 ), which was involved in response to hypoxia (GO:0001666) (Supplementary Table S6).For low-altitude cattle, the heart-related module (MElightyellow, r = 0.99, p = 1 × 10 -15 ) was enriched in cardiac right ventricle morphogenesis (GO:0003215), negative regulation of cardiac muscle cell apoptotic process (GO:0010667), and heart development (GO: 0007507).Lung-related module (MEdarkorange, r = 0.99, p = 7 × 10 -16 ) was enriched in lung development (GO:0030324).Additionally, in the tissue-related module of each species, the terms related to bile secretion were found in the liver-specific module, functions related to ion transport and absorption are enriched in the kidney-specific module, and immune and skeletal muscle cell differentiation-related functions were enriched in the spleen-specific and muscle-specific module, respectively (Supplementary Table S6).
Furthermore, genes with the highest intramodular connectivity were selected as hub genes in each tissue-related module in both gayal and yak.We identified 4-51, 6-44, and 2-46 hub genes in six tissue-related modules of gayal, yak, and cattle, respectively (Supplementary Table S7).Among the heart-related modules of the three species, there are fourteen common hub genes between gayal and yak, of which five genes are involved in cardiac contraction (bta04260), including UQCRC1, UQCRQ, UQCRFS1, UQCR11, and COX5A.Among them, UQCRC1, UQCRQ, UQCRFS1, and UQCR11 are important subunit of mitochondrial complex Ⅲ (Fernández-Vizarra and Zeviani, 2015).Mitochondria play a key role in cardioprotection and heart contraction energy needed for pumping blood to oxygenate the body organs (Lemieux and Hoppel, 2009).UQCRC1 directly affects mitochondria function related to cardioprotection (Yi et al., 2020).COX5A is the component of the cytochrome c oxidase, which is central to oxidative phosphorylation and ATP generation (Huang et al., 2019).However, no common hub gene was found between gayal and low-altitude cattle, and only four hub genes are shared by yak and cattle, including KCNJ8, MTUS2, P4HTM, and BTG3.KCNJ8 encodes a subunit of an ATPsensitive potassium channel and regulates potassium inward rectifier, which is related to heart development (Erginel-Unaltuna et al., 1998).There are one (COX6A2) and two (CYC1 and UQCR10) non-common hub genes in the heartrelated modules of gayal and yak, respectively, which are involved in cardiac contraction.In addition, four hub genes (HAND2, CHRM2, EDNRA, and ZFPM2) associated with heart development and regulation were uniquely identified in the greenyellow module of the yak heart.Among these, CHRM2 was an up-regulated DEG in the yak heart (log 2 fold changes of 1.14), which plays a fundamental role in autonomic regulation of the heart, such as heart rate recovery after maximal exercise (Hautala et al., 2006).HAND2 is essential for cardiac morphogenesis, vascular development, and the regulation of angiogenesis (Anderson et al., 2016).EDNRA encodes the endothelin 1 (EDN1) receptor, a gene involved in the vasoconstriction mechanism (Sharma et al., 2014), as well as closely related to pulmonary hypertension and HIF activity, which is positively selected in Tibetan (Simonson et al., 2010).ZFPM2, regulates the activity of GATA family proteins, plays a vital function in cardiac morphogenesis and coronary vascular development (Luo et al., 2020).
Among the hub genes of the lung-related module, six common hub genes were identified among gayal and yak, including CAPS, TSNAXIP1, GPRIN2, SPEF1, BMP3, and EFHC1.Gayal and yak have two (BMP3 and SEC14L3) and three (BMP3, DNAAF1, and N4BP3) hub genes in common with low-altitude cattle respectively.Among them, BMP3, a hub gene shared by all three bovine subfamily species, has been reported to play a regulatory role in the morphogenesis and/or function of the human lung (Vukicevic et al., 1994).CAPS encodes a calciumbinding protein and plays a role in the regulation of Ca 2+ transport (Pan et al., 2019), which may be correlated with pulmonary artery smooth muscle contraction (Waypa and Schumacker, 2002).SEC14L3 is an up-regulated DEG in the gayal lung and plays an important role in maintaining the homeostasis of airway epithelial cells (Shan et al., 2009).DNAAF1 gene is associated with lung development (Duquesnoy et al., 2009).As for the non-common hub genes in the lung-related module of gayal, we discovered several genes (EDN3, GJA5, and CCDC39) related to the tracheal epithelium and pulmonary vasoconstriction.Such as EDN3, an endotheliumderived vasoconstrictor peptide, which was up-regulated in the lungs of gayal, and may play an important role in enhancing pulmonary vasoconstrictor response (Masaki, 2000).GJA5 is the predominant gap junction protein present in vascular endothelium, plays an important role in coupling between cells in the vascular wall (Ebong et al., 2006).CCDC39 was related to lung development (Blanchon et al., 2012).Therefore, these genes related to promoting cardiac blood circulation and pulmonary vasoconstriction may play crucial roles in the environmental adaptability of gayal and yak.

DISCUSSION
In this study, we explored the high-altitude adaptability of Bovidae using the transcriptome data of six tissues from gayal and published data on yak and lowland cattle (Tang et al., 2017).Compared with numerous previous studies that only used the same reference genomes (yak or cattle genome) for different bovine subfamily species (Tang et al., 2017;Lan et al., 2018;Xin et al., 2019;Wang et al., 2020), here we optimized the experimental design, using the corresponding reference genomes for three species to reduce the gene sequence differences of different species due to phylogenetic evolution, and only retained reads that mapped to the 1:1 orthologous gene which had high-confidence alignments across the three genomes for subsequent analysis.To reduce the difference in the sequencing depth of transcripts from different sequencing platforms, we used the TMM standardized method of edgeR for correction with reference to the method of Blake et al. (Blake et al., 2020).Consistent with previous studies (Necsulea and Kaessmann, 2014;Tang et al., 2017;Hao et al., 2019), we found that the differences in global gene expression among tissues were more significant than that among altitudes or the same tissue of different species.Among them, four tissues (lung, liver, spleen, and kidney) showed tissue-specific expression patterns, while bovine species-specific expression patterns were shown in the tissue of heart and muscle (Necsulea and Kaessmann, 2014;Tang et al., 2017;Hao et al., 2019).In addition, we identified numerous DEGs in gayal and yak compared with low-altitude cattle, which might have an important role in high altitude adaptation.WGCNA was further used to explore the core genes of each tissue that regulate tissue function in the high-altitude adaptation.

Expression Regulation of Genes Involved in the Regulation of Blood Cell Development
Hypoxia is one of the most serious challenges faced by highaltitude animals, and oxygen supply is directly related to the development of red blood cells (Zhang D. et al., 2017).Studies in mice (Li et al., 2011) and yak (Xin et al., 2019) have shown an increase in the number of red blood cells and platelets under hypoxia, but the number of granulocyte-macrophage progenitor cells periodically declined.In this study, similar results were found in the gayal.The hematopoietic cell lineage pathway was uniquely enriched in multiple tissues (heart, kidney, lung, and spleen) of gayal vs. the cattle group.Among these, IL11RA is a member of the hematopoietic cytokine receptor family (Ng et al., 2021), which was uniquely up-regulated by 1.33 times in the heart of gayal.GP1BA is the alpha subunit of platelet surface membrane glycoprotein Ib (Luo et al., 2007), were 9.58, 3.90, 7.77, 11.41, 12.38, 12.39-fold higher in the tissues (corresponding to heart, spleen, lung, kidney) of gayal than in the cattle, suggesting that the number of red blood cells and platelets might be increased.As for the regulatory mechanisms underlying granulocytemacrophage progenitor cells, we observed three genes (IL1B, CSF1, and IL6R) were down-regulated in multi tissues of gayal, which might be a reason for the decreased proportion of granulocytes.Among them, IL1B is produced by activated macrophages as a proprotein (Lappas, 2013), which was downregulated 1.94-fold in the lung of gayal.CSF1 controls the production, differentiation, and function of macrophages (Hamilton et al., 2016), displaying log 2 fold changes of −1.06, −1.10, and −1.76 in the heart, spleen, and muscle, respectively.In addition, Interleukin 6 is involved in the maturation of B cells.IL6R encodes a subunit of the interleukin 6 (IL6) receptor complex (Liu et al., 2017) and was found down-regulated in all tissues of gayal, which suggests lymphocytes may decrease, which is consistent with the previous results on the hypoxia adaptation of yak (Xin et al., 2019).
In addition, in the high-altitude related gene set that we searched (Supplementary Table S5), several crucial positive selection genes related to erythrocyte development were found differential expression in the gayal and yak.EGLN1, a key positive selection gene in the HIF pathway (Bigham et al., 2010), was found consistently down-regulated in all tissues of gayal, which contribute to enhancing HIF1-α activity and has a role in increased Hb levels via erythropoiesis under hypoxia conditions (Epstein et al., 2001;Xiang et al., 2013).HFE is a positively selected gene in Tibetans (Yi et al., 2010), which indirectly regulates iron balance by reducing the entry of transferrin into cells (Salter-Cid et al., 1999).Compared with lowland cattle, the appropriate down-regulation of HFE gene expression in gayal can accelerate iron storage to facilitate the synthesis of hemoglobin to cope with the reduced oxygen partial pressure (Muckenthaler et al., 2020).In the high-altitude related DEG set of yak, we also found EGLN1 was consistently downregulated in all tissues of yak.In addition, HOXB6 is considered to be a marker for erythropoiesis (Magli et al., 1997), and its expression was uniquely up-regulated 1.85 and 1.11-fold in the spleen and kidney of the yak, respectively, which may contribute to yak adaptation to the hypoxia environment (Xin et al., 2019).

Expression Regulation of Genes Involved in Angiogenesis
The regulation of angiogenesis by hypoxia links blood oxygen supply to metabolism (Pugh and Ratcliffe, 2003).In this study, the shared and unique DEGs involved in blood vessel and vascular development were identified.For example, glutamyl aminopeptidase (ENPEP), a member of the M1 family of endopeptidases involved in blood pressure regulation and blood vessel formation (Holmes et al., 2017), was up-regulated in the lungs of gayal and yak with the log 2 fold changes of 2.78 and 2.98, respectively.Additionally, we found that SRPX2 and ESM1 were up-regulated in the spleen of gayal and yak relative to lowland cattle, and SAT1 and SERPINE1 genes were up-regulated in the kidney and muscle of gayal and yak relative to lowland cattle, respectively, which are angiogenesis promoting genes (Aitkenhead et al., 2002;Mohr et al., 2017;Xiao et al., 2020).HEY1 is the downstream effector of Notch signaling required for cardiovascular development (Aitkenhead et al., 2002) with log 2 fold changes of 1.08 and 1.71 in the heart and lung of gayal, respectively, while displaying log 2 fold changes of 1.25 and 1.03 in the yak, respectively.In addition, MEOX2 has been reported to have a role in inhibiting angiogenesis (Dhahri et al., 2020), and was down-regulated in all 6 tissues of gayal and the kidney of yak.Research has found that VHL-deficient mice show increased vasculogenesis in embryonic cells (Wang et al., 2017), and the VHL gene was down-regulated in expression in all six tissues of yak and down-regulated in the lung of gayal.In addition, we have identified two unique DEGs (RAMP2 and TGFBI) in gayal that promote angiogenesis.RAMP2 mediates the pro-angiogenic effect (Albertin et al., 2010), which was up-regulated in all six tissues of gayal.TGFBI is a matricellular protein-coding gene that plays an important role in tumor angiogenesis (Fico and Santamaria-Martínez, 2020), which was up-regulated in the liver and kidney of gayal with log 2 fold changes of 1.34 and 1.45, respectively.Additionally, DEGs (THY1 and PLCD3) that promote angiogenesis were also uniquely found in yak.THY1 was expressed on mouse tumor-associated lymphatic vessels and blood vessels (Jurisic et al., 2010), displaying log 2 fold changes of 1.63 and 1.03 in the spleen and lung of yak.PLCD3 induce angiogenesis in human endothelium (Kim et al., 2002), which was up-regulated in all six tissues, displaying log 2 fold changes of 1. 77, 1.83, 2.07, 3.20, 2.12, and 1.15 in the heart, liver, spleen, lung, kidney, and muscle of yak, respectively.This suggests that angiogenesis may be conducive to the adaptation of gayal and yak to a high-altitude environment.

Expression Regulation of Genes Involved in the Immune System
Increasing evidence shows that the harsh environment of hypoxia may regulate the immune system (Mishra and Ganju, 2010;Gaur et al., 2020).In the present study, we found that the immune system of the spleen in gayal and yak was highly regulated.DEGs of the spleen in gayal were significantly enriched in lymphocytic chemotaxis, inflammatory response, and cytokine receptor interaction pathways (Supplementary Table S4), and the DEGs of the spleen in yak were also significantly enriched in TNF signaling pathway, cytokine receptor interaction, and chemokine signaling pathway (Supplementary Table S4).In the spleen of gayal vs. cattle group, chemokines and their receptors like CCL4, CCL5, CCL14, CCL16, and CXCR4 displayed log 2 fold changes of 2.33, 2.78, 1.24, 1.83, and 1.18, respectively, which mediates chemokinetic inflammatory response (Loetscher et al., 1994;Ariel et al., 2000); Interleukins like IL12RB2 promotes the proliferation of T-cells as well as NK cells (Presky et al., 1996), IL23A contribute to the production of proinflammatory cytokines (Hoeve et al., 2006), and IL12RB2 and IL23A with log 2 fold changes of 1.77 and 1.09, respectively; tumor necrosis factors like TNFRSF8, TNFRSF4, and LTBR displayed log 2 fold changes of 4.06, 1.41, and 1.32, respectively, which activated T and B cell expression (Nishikori et al., 2005).CD40LG regulates B cell function by engaging CD40 on the B cell surface (Henn et al., 1998), which was up-regulated 2.63-fold in the spleen of gayal.In the spleen of the yak vs. cattle group, CXCR4 and CD40LG were also up-regulated, with log 2 fold changes of 1.52 and 2.83, respectively.CXCL14 which fulfills a unique role in antimicrobial immunity was up-regulated 2.15-fold in the spleen of yak (Maerki et al., 2009).These changes might help gayal and yak resist inflammation and disease in high-altitude harsh environments.
High altitude pulmonary edema (HAPE) is a serious lifethreatening disease in humans and animals characterized by uneven vasoconstriction of pulmonary blood vessels due to hypoxia (Stream and Grissom, 2008).It has been confirmed that HAPE is related to coagulation activation and fibrin formation enhancement (Mannucci et al., 2002).In addition, hypoxia can lead to increased platelet activity, significantly increasing the risk of thrombosis (Shilei et al., 2016).In this study, we found that a large number of DEGs are related to platelet activation and blood coagulation in the two comparison groups.Coagulation factors play an important role in the endogenous and exogenous blood coagulation process in animals (Thiel et al., 2017).In this study, compared with lowaltitude cattle, the expression levels of coagulation factor II (F2) and coagulation factor V (F5) in the lung of gayal were separately reduced 1.26 and 1.64-fold, and displayed log 2 fold changes of −2.33 and −1.23 in the lung of yak, respectively.In addition, coagulation factor X (F10) was found uniquely down-regulated 5.62 times in the lungs of yak.Meanwhile, serine protease inhibitors are inhibitors of blood coagulation factors (Seixas et al., 2007), and we have found the expression level of SERPINA5, SERPINF2 and SERPING1 were uniquely increased in the lung of gayal, and the corresponding log 2 fold changes of 1.80, 3.06, and 1.48, which may play an important role in the mechanism of pulmonary edema resistance, thus reducing the risk of pulmonary edema.In addition, impaired decomposition of bradykinin and its metabolites may be an important cause of pulmonary edema (De Maat et al., 2020).We found that the expression of the kininogen 1 (KNG1) gene was down-regulated by 2.46 times in the lung of yak, while down-regulated by 12.39 and 7.38-fold in the heart of gayal and yak than that of lowland cattle, respectively.These results suggest yak and gayal have developed a mechanism to prevent the occurrence of severe pulmonary edema.

Expression Regulation of Genes Involved in Energy Metabolism
High energy metabolism or reduced energy turnover is essential for survival under hypoxia (Qiu et al., 2012).In this study, mitochondrion-related categories were widely enriched in the tissues of gayal and yak (Supplementary Table S4).Mitochondria is the main place for aerobic respiration of cells, which provides ATP for organisms through the oxidation of sugar, fats, and amino acids (Solaini et al., 2010).Similarly, several pathways related to carbohydrate and fat energy metabolism were enriched, such as glycolysis/gluconeogenesis in yak lung, oxidative phosphorylation, and galactose metabolism in multiple tissues of gayal, and PPAR signaling pathway in the heart of gayal and yak.Among the genes related to these functions, TPI1 is essential for efficient energy production and plays an important role in glycolysis (Shimoda et al., 2012).G6PC (glucose-6-phosphatase) is a key enzyme of gluconeogenesis (Jia et al., 2012).We found that the gene expression of TPI1 was upregulated in all tissues of both gayal and yak, and G6PC was down-regulated in three tissues (heart, spleen, and lung) of both gayal and yak, which is conducive to the production of ATP through strengthened glycolysis and reduced gluconeogenesis.Oxidative phosphorylation is an important process for ATP syntesis in aerobic organisms.We found that most genes related to oxidative phosphorylation were up-regulated in gayal, such as COX6A2, NDUFA9, COX8B, NDUFB7, NDUFA11, NDUFB5, UQCR10, NDUFS8, UQCRQ, PPA1, and NDUFS3.In addition, the significant down-regulation of COX5B found in all tissues of gayal and yak may be beneficial to hypoxia adaptation (Trueblood et al., 1988).Additionally, galactose metabolism (bta00052) was uniquely enriched in the tissues (lung, liver, spleen, and muscle) of gayal.Galactose is a kind of monosaccharide, most of which are converted into glucose in the liver, and then incorporated into glycogen or used for energy metabolism by glycolysis.Lactase (LCT), a key enzyme involved in the conversion of galactose to glucose was up-regulated in all six tissues of gayal.The PPAR signaling pathway, a classic pathway in lipid catabolism was significantly enriched in the heart of gayal and yak.EHHADH, a dehydrogenase in the βoxidase system (Yeh et al., 2006), was down-regulated 4.02 and 2.69 times in the heart of gayal and yaks.Additionally, the utilization of fatty acids first needs to be activated to produce fatty acyl CoA under the catalysis of the ACSL and SLC27A family (Bowman et al., 2016).In this study, we found that SLC27A2 and SLC27A6 of the SLC27A gene family were down-regulated 6.24 and 3.24 times in the heart of gayal, and SLC27A2 and SLC27A4 were down-regulated 6.81 and 2.08 times in yak heart, respectively.And the ACSL3 and ACSL4 genes expression of ACSL gene family members were found uniquely decreased in the heart of gayal.Therefore, our results suggest gayal and yak may adapt to high altitude environments by strengthening carbohydrate metabolism and reducing fatty acid degradation although the related regulation genes may be different between gayal and yak.
Despite the fact that we found numerous genes that may be related to the high-altitude adaptation of gayal and yak, we must admit that the limitation of age differences between gayal and other species, the gayal being a younger species, may affect the experimental results obtained to some extent.Therefore, it is necessary to conduct comparative transcriptome analysis on more bovine species of the same age that are combined with physiologic experiments to verify our findings, and better clarify the high-altitude adaptation mechanism of bovine subfamily species.

CONCLUSION
This study shows a comprehensive multi tissues transcriptome expression landscape of high-and low-altitude bovine populations and reveals the gene expression regulation associated with highaltitude acclimatization.The gene expression profiles of gayal, yak, and cattle showed tissue-specific expression patterns.The comparative analysis of six hypoxia-related tissues of high-and low-altitude bovine subfamily species highlights numerous DEGs and terms underlying associated with hypoxia.Notably, categories related to angiogenesis, blood coagulation, energy metabolism, and the immune system were commonly enriched in the tissue DEGs of gayal vs. cattle group and yak vs. cattle group, and we found that many expression regulatory genes related to these functions in gayal and yak are different, which may serve as an important regulatory mechanism for gayal and yak to adapt to the corresponding local environment.In addition, we also found numerous hub genes related to myocardial contraction, energy metabolism, and pulmonary vasoconstriction by WGCNA analysis.Overall, our study lays a foundation for further study on the environmental adaptability of mammals in subtropical plateau and provides a valuable basis for resource utilization in high-altitude areas.

FIGURE 1 |
FIGURE 1 | Gene expression patterns across the three bovine species: yak, gayal, and cattle.(A) symmetrical heatmap of Pearson's correlation coefficients between all pairs of samples across all genes.The key is hierarchically clustered (complete method, correlation distance).(B) PCA of the log-transformed normalized expression levels of all orthologs across all species and tissues.Species are represented by point shape; tissues are represented by point color.

FIGURE 2 |
FIGURE 2 | DEGs upset map of all DEGs in six tissues of (A) gayal vs. cattle group and (B) yak vs. cattle group.

FIGURE 3 |
FIGURE 3 | GO functional and KEGG pathway enrichment analysis of DEGs in the (A) heart, (B) liver, (C) spleen, (D) lung, (E) kidney, and (F) muscle for the gayal vs. cattle group.

FIGURE 4 |
FIGURE 4 | GO functional and KEGG pathway enrichment analysis of DEGs in the (A) heart, (B) liver, (C) spleen, (D) lung, (E) kidney, and (F) muscle for the yak vs. cattle group.

FIGURE 5 |
FIGURE 5 | Venn diagram indicating differentially expressed genes that were shared among the two high-and low-altitude pairs in (A) heart, (B) liver, (C) spleen, (D) lung, (E) kidney, and (F) muscle.Numbers in red and blue indicate genes up-and down-regulated in the high-altitude Bovina species relative to low-altitude cattle.

FIGURE 6 |
FIGURE 6 | WGCNA analysis of gayal tissue samples.(A) Hcluster diagram of gayal tissue samples.(B) Plot of scale-free topology and mean connectivity in regard to soft-thresholding power for samples.(C) Hierarchical clustering tree (dendrogram) of genes based on co-expression network analysis.(D) Heatmap of the correlation between the module eigengenes and tissues of gayal.

FIGURE 7 |
FIGURE 7 | WGCNA analysis of yak tissue samples.(A) Hcluster diagram of yak tissue samples.(B) Plot of scale-free topology and mean connectivity in regard to soft-thresholding power for samples.(C) Hierarchical clustering tree (dendrogram) of genes based on co-expression network analysis.(D) Heatmap of the correlation between the module eigengenes and tissues of yak.

FIGURE 8 |
FIGURE 8 | WGCNA analysis of cattle tissue samples.(A) Hcluster diagram of cattle tissue samples.(B) Plot of scale-free topology and mean connectivity in regard to soft-thresholding power for samples.(C) Hierarchical clustering tree (dendrogram) of genes based on co-expression network analysis.(D) Heatmap of the correlation between the module eigengenes and tissues of cattle.

TABLE 1 |
The number of DEGs identified in six tissues for gayal and yak compared with low-altitude cattle.