Six tissues (heart, kidney, liver, lung, skeletal muscle, and spleen) were collected from three adult female gayal living in the semi-wild 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.

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 first-strand complementary DNA (cDNA). Subsequently, the second-strand 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.

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.

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.

After performing normalization, we obtained the TMM-normalized log₂-transformed CPM values for each gene. The principal component analysis (PCA) was performed with the prcomp function in R using normalized log₂-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.

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₂ 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).

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).

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.55 G 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.

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).

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 (species-specific 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.

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 up-regulation 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 six tissues of gayal and yak compared with lowland cattle, respectively (Figure 5). By functional annotation analysis, we observed the shared DEGs of the heart were enriched in the oxidation-reduction process (GO:0055114), and positive regulation of blood coagulation (GO:0030194). Lung-shared DEGs were involved in platelet activation (GO:0030168), and 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).

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).

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 ² 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 tissue-related module.

Indeed, the relevant modules enriched in six tissues are highly related to the tissue-specific functions. For example, in the tissue-related 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 ATP-sensitive 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 heart-related 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₂ 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 calcium-binding protein and plays a role in the regulation of Ca²⁺ 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 endothelium-derived 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.

Hypoxia is one of the most serious challenges faced by high-altitude 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 granulocyte-macrophage 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 down-regulated 1.94-fold in the lung of gayal. CSF1 controls the production, differentiation, and function of macrophages (Hamilton et al., 2016), displaying log₂ 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 down-regulated 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).

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₂ 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₂ fold changes of 1.08 and 1.71 in the heart and lung of gayal, respectively, while displaying log₂ 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₂ 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₂ 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₂ 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.

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₂ 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₂ fold changes of 1.77 and 1.09, respectively; tumor necrosis factors like TNFRSF8, TNFRSF4, and LTBR displayed log₂ 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₂ 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 life-threatening 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 low-altitude 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₂ 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₂ 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.

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 up-regulated 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.