Whole-transcriptome analysis of heart in bovidae reveals regulatory pathways associated with high-altitude adaptation

Background: the yak (Bos grunniens) is one of the major livestock that can survive the extremely cold, harsh and oxygen-poor conditions and provide meat, milk and transportation for Tibetans living in the Qinghai-Tibetan Plateau, and there are some cattle migrate to yak habitat but cannot survive the extremely conditions well. However, the regulatory mechanisms that drive hypoxic adaptation of yak remain elusive. Results: Thus, we collected heart tissues from Leiwuqi (LWQ) yak and their relatives, the migration cattle in LWQ yak habitat (HAC) and low-altitude cattle (LAC), respectively, for RNA sequencing, 178 co-differentially expressed protein-coding transcripts (CETs) including 6 down-regulated and 172 up-regulated CETs were discovered in both LAC-vs-LWQY and LAC-vs-HAC comparison groups. There were 133, 11 and 3 lncRNA transcripts were specific expressed in LWQY, HAC, and LAC groups, respectively. 2 and 230 lncRNA transcripts differentially expressed in LAC-vs-LWQY and LAC-vs-HACrespectively, but no lncRNA transcripts co-differentially expressed in those two comparison groups. Among the 58 miRNAs that were co-differentially expressed, 18 were upregulated and 40 were down-regulated. In addition, 640 (501 up-regulated and 139 down-regulated) and 152 (152 up-regulated and 1 down-regulated) circRNAs showed differential expression in LAC-vs-LWQY and LAC-vs-HAC comparison groups, respectively, and 53 up-regulated co-differentially expressed circRNAs were shared. Multiply CETs which are the target genes for miRNAs/lncRNAs were significantly enriched in high-altitude adaptation related processes, for instance, T cell receptor signaling, VEGF signaling, and cAMP signaling. A competing endogenous RNA (ceRNA) network by integrating competing relationships among CETs, miRNAs, lncRNAs and circRNAs were constructed. Furthermore, we constructed a hypoxic adaptation related ceRNA network to identified that 8 miRNAs and 15 circRNAs were predicted to regulated MAPKAPK3, PXN, NFATC2, All these results indicated that and sever environment in yak. A total of 34218 circRNAs were obtained in the nine libraries using the CIRI software, 15283 (66.38%) of them were generated from exons of protein-coding genes, 1102 (4.79%) were from intergenic regions, 1813 (7.87%) were from exon-intron, and 1739 (7.55%) were from introns, the remaining (3077, 13.36%) were antisense, and the average length of circRNAs was 2404 bp in cattle. While, 12293 (62.97%) of them were generated from exons of protein-coding genes, 1610 (8.25%) were from intergenic regions, 2308 (11.82%) were from exon-intron, and 854 (4.37%) were from introns, the remaining (2475, 12.68%) were antisense, and the average length of circRNAs was 2719 bp in yak (Supplementary file16). In addition, 19522, 25943, and 11137 circRNAs were obtained in LWQY, HAC and LAC group, respectively (Supplementary file16). These results suggested that the circRNAs in bovidae hearts were from diverse genomic regions, and the average length in different species was different. Comparing the whole transcriptomes of yak and cattle can provide insight into the regulatory mechanisms of high-altitude adaptation at a comprehensive view, including among mRNAs, miRNAs, lncRNAs, and circRNAs. Here, we performed whole transcriptomes analysis in the heart of Leiwuqi (LWQ) yak and their relatives, the migration cattle in LWQ yak habitat (HAC) and low-altitude cattle (LAC), respectively, and identified comparable numbers of transcriptomes species using yak and cattle genome, respectively. In addition, we also constructed the hypoxia-adaptation ceRNA network. We identified 178 differentially expressed genes DEGs between LAC-vs-LWQY and LAC-vs-HAC group, further examination revealed that 172 DEGs are up-regulated in LAC-vs-HAC groups. KEGG analysis reveals these DEGs were enriched in hypoxia-adaptation related pathways, in which some important oxygen transportation-, environmental sense-, energy metabolism-, and immune-relevant genes, such as COX7C , MAPKAPK3 , PIK3R5 , NFATC1 and NFATC1 were significantly up-regulated in HAC and LWQY compared to LAC. COX7C is one subunits of cytochrome c oxidase (COX), which is the terminal component of the mitochondrial respiratory chain. As a regulator in hydrogen ion transmembrane transport, COX7C play an important role in energy metabolism. MAPKAPK3 is a serine/threonine protein kinase of the p38 signaling pathway that is activated by a variety of stress stimuli and is implicated in cellular responses and gene regulation during the adaptation to hypoxia. PIK3R5

The yak is an indigenous and rare bovine species that distributed on the Qinghai-Tibet Plateau and adjacent areas at altitudes above 2500m. Archaeological evidence indicates that yak and cattle, both originated from wild origin cattle, are two closely related species which diverged from five million years ago [1][2]. Previous studies show that the genomes of these two species have strong similarities; including both of them have 30 chromosomes and similar karyotypes [3], and the extensive synteny of the two ruminants' genomes [2]. Yak but cattle suffer from sever pulmonary hypertension when migrated to the yak habitat [4][5]. After evolved over million years, yak have developed many anatomical and physiological traits that enable them to survive at high altitudes, for instance, larger lungs and hearts [3], shorter tongue [6], stronger environmental sense [3], higher energy metabolism [3,7], and lack of hypoxic pulmonary vasoconstriction [3]. Hence, it is necessary not only to use the migration and low-altitude cattle to study the hypoxic adaptation mechanism, but also to use yak. Recent studies on morphology [8], anatomy [9], physiology [10], genomics [3], mRNA [11][12] and small RNAs [13] of yak have focused on the high-altitude adaptation mechanism, but limited researches are centered on the whole transcriptome and their regulatory network in high altitude adaptation.
The universal whole transcriptome refers to the clustering of all transcripts in tissue or cells, including mRNA and non-coding RNA (ncRNA) [14]. Nowadays, among the different kinds of ncRNAs, microRNA (miRNA), long noncoding RNA (lncRNA) and circular RNA (circRNA) receives researcher's most attention, as all of whose regulatory functions are associated with mRNA [15,16]. MiRNA and lncRNA are two classes of endogenous, linear ncRNAs, which play important roles in various biological and pathological processes, for instance, hypoxia adaptive response [17,18], lipid and glucose metabolism [19,20], carcinogenesis and tumor progression [21]. However, the regulatory mechanism of these two kinds of ncRNAs is different. MiRNAs generally post-transcriptionally down regulate mRNA stability or block their translation by binding to the 3'-untranslated region (3'-UTR) of specific mRNA targets in plants and animals [22]. It has been proved that there are 70 miRNAs differentially expressed in lung between yak and cattle, and those miRNAs are associated with hypoxia-related functions [13]. lncRNA, as a transcription regulator, can not only through chromatin modification, but also via base pairing of mRNA and using decoys of RNA-binding proteins/miRNAs to affect gene expression [23]. Previous studies shows that histone methylases MLL1 coordinates with hypoxiainducible factors 1α(HIF1α) and histone acetyltransferase p300 and regulate hypoxia-induced lncRNA HOTAIR expression [24]. For circRNA, recent studies show that it can inhibit miRNA functions as miRNA sponges, and modulate gene expression both transcriptionally and posttranscriptionally (such as mRNA splicing and stability) [25]. Furthermore, circRNAs can interact with proteins to regulate gene expression. To explain the correlation of genome size and organism complexity, Salmena et al. raised "ceRNA hypothesis" which hypothesize that mRNA, transcribed pseudogenes and lncRNA "talk" to each other through miRNA response elements (MREs), then forms large-scale regulatory networks across the transcriptome [26]. Nowadays, numerous studies have verified these regulatory networks really exist and play pivotal roles in almost all the pathological and biological processes.
LWQ yak are distributed in Changdu city (Tibet province, China) with the altitude above 3600m. In the present study, we used heart (the hypoxia-sensitive tissues) tissues from LWQ yak, the migration cattle in LWQ yak habitat and low-altitude cattle, to construct RNA library respectively, and performed high throughput sequencing to investigate: (Ⅰ) whether ncRNAs were involved in the mechanism of high altitude adaptation; (Ⅱ) the high altitude adaptation mRNAs and ncRNAs; (Ⅲ) which pathways those differentially expressed ncRNAs enriched; (Ⅳ) whether those ncRNAs "talk" to each other.

Overview of RNA-sequencing
To obtain a global view of the yak and cattle heart transcriptome (including mRNA, lncRNA, miRNA, and circRNA) and identify whether they are related to high altitude adaptability. Two kinds of cDNA libraries, a strand-specific library and a miRNA library, were constructed and sequenced for each sample. A total of 973,288,574 raw reads were obtained from strand-specific libraries, with an average of 108,143,174 in each library. After a series of strict criteria, a total of 959,028,938 raw reads were defined as high quality reads for identifying the mRNAs, lncRNAs and circRNAs  detected between LAC and HAC libraries, and 4,742 events were significant between these two groups. The distribution of alternative splicing events in the two groups is similar, with three key type events being prevalent including skipped exon (SE), alternative 3' splice site (A3SS) and mutually exclusive exon (MXE) (Fig1A), implying that vast novel protein-coding transcripts identified were mainly generated by these alternative splicing events, which play critical roles in the hypoxic adaptation and high-altitude diseases of cattle.  Fig 2) were carried out to identify biological functions of the 178 DETs with the cattle reference genome. Among these, the most important cellular components involved the cell, cell part, organelle, and extracellular matrix component. The molecular functions consisted of binding, catalytic activity, nucleic acid binding transcription factor activity, and signal transducer activity. In biology processes, the important GO terms were biology adhesion, immune system process, response to stimulus, and cell adhesion. In the KEGG pathway analysis, numerous genes were enriched from signal transduction pathways including the immune system, and signal transduction system. In the enrichment of signal pathways, we identified VEGF signaling pathway [27], cAMP signaling pathway [28], oxidative phosphorylation, B/T cell receptor signaling pathway, natural killer cell mediated cytotoxicity pathway [29], and cardiac muscle contraction [30] were related to hypoxic adaptation system that is important for oxygen transportation, environmental sense, energy metabolism, and infection and immune defense. There were 15 genes were associated with the above pathways, including cytochrome c oxidase subunit 7C (COX7C), mitogen-activated protein kinase-activated protein kinase 3 (MAPKAPK3), BCL2 associated agonist of cell death (BAD), vav guanine nucleotide exchange factor 1 (VAV1), nuclear factor of activated T cells 2 (NFATC2), nuclear factor of activated T cells 1 (NFATC1), and phosphoinositide-3-kinase regulatory subunit 5 (PIK3R5), and these genes owned significant higher (> 10-fold) expression level in the heart tissues of HAC and LWQY compared to that of in LAC. These results indicated that the mechanism of high-altitude adaptation involve diverse complex cellular processes and changes in related genes and kinases, which are worthy of future research efforts.
We also detected involvement of a few interesting pathways involved in endocrine system, and metabolism of other amino acids, such as the thyroid hormone signaling pathway [31], D-Glutamine and D-glutamate metabolism [32], and previously studies shown that these pathways were associated In LWQY group, the average expression of lncRNA transcripts was lower than that of protein-coding transcripts in the three groups (Fig 3A), the average length and number of exons of protein-coding transcripts were longer than the lncRNA transcripts ( Fig 3C, 3E), and the number of exons of lncRNAs and proteincoding transcripts were concentrated at 2 and 5 exons, respectively ( Fig 3E). The similar results were observed in cattle ( Fig 3B, 3D,3F). In addition, the average GC content of protein-coding and lncRNA transcripts in LWQY was 51.12% and 49.41% (Fig 3G), respectively. While in cattle the GC content were 48.79% and 51.38%, respectively ( Fig 3D). The present results that are in accordance with those of previous studies [33][34][35].
In total, there were 2 and 230 lncRNA transcripts differentially expressed in LAC-vs-LWQY and LAC-vs-HAC group (Fig 2E), but no lncRNA transcripts co-differentially expressed in those two comparison groups with the limited differentially expressed lncRNAs in LAC-vs-HAC comparison group. Through cis-, trans-, and antisense-regulatory relationship analysis, we detected 168 potential target proteincoding genes of the differential lncRNA transcripts. Functional analysis showed that these target  (Table 2). A total of 299 and 90 miRNAs were differentially expressed in LAC-vs-LWQY and LAC-vs-HAC groups, respectively (Supplementary file11 and S12). 58 miRNAs were shared co-differentially expression in these two comparison groups, including 18 were up-regulated and 40 were down-regulated ( Fig 1D).
To further investigate the functions of these differentially co-expressed miRNAs, their target genes were predicted and function annotation were performed. A total of 43930 target protein-coding genes for co-differentially expressed miRNAs were identified (Supplementary file13). Among them, 125 genes had transcripts identified as CETs and also significantly negatively correlated with miRNA expression (Supplementary file14), and were assigned as intersection genes, which were more likely to be predicted as miRNA target genes in heart of Bovidae. Similar to the mRNA expression results, the largest numbers of target genes were clustered in biological processes, followed by molecular function and cellular component (Supplementary file15). KEGG analysis revealed that these intersection genes were involved T cell receptor signaling, VEGF signaling, and cAMP signaling (Fig 4).
All these results indicated that miRNAs regulated a variety of processes to adapt the hypoxic, cold and sever environment in yak.
The circRNAs expression analysis and circRNA-miRNApathway relationship network construction A total of 34218 circRNAs were obtained in the nine libraries using the CIRI software, 15283 (66.38%) of them were generated from exons of protein-coding genes, 1102 (4.79%) were from intergenic regions, 1813 (7.87%) were from exon-intron, and 1739 (7.55%) were from introns, the remaining (3077, 13.36%) were antisense, and the average length of circRNAs was 2404 bp in cattle. While, 12293 (62.97%) of them were generated from exons of protein-coding genes, 1610 (8.25%) were from intergenic regions, 2308 (11.82%) were from exon-intron, and 854 (4.37%) were from introns, the remaining (2475, 12.68%) were antisense, and the average length of circRNAs was 2719 bp in yak were shared (Fig 1E). It has been reported that circRNAs can bind to miRNAs and consequently repress their function. To determine whether the 53 differentially expressed cirRNAs could affect gene post-transcriptional regulation by binding to miRNAs, we use the 3 software to predict their target miRNAs. Interesting, these 53 differentially expressed cirRNAs could bind to 881 miRNAs, and 57 of them had transcripts identified as co-differentially expressed miRNAs and also significantly negatively correlated with 125 DETs (Supplementary file19), which suggested that cirRNAs may play important roles as miRNA sponges across a wide array of biological processes in the heart to protect yak from the harsh environment of the plateau.

Construction of ceRNA network
Based on the data for co-differentially expressed protein-coding transcripts, circRNA transcripts and miRNA transcripts in LAC-vs-LWQY and LAC-vs-HAC groups, we used the three softwares, mireap, miRanda and TargetScan, to identify biological targets of each miRNA from the protein-coding and circRNA transcripts that showed a significantly negative correlation with miRNA expression, subsequently obtained the mRNA-miRNA and circRNA-miRNA pairs, and reserved the mRNA-miRNA pairs that the mRNA were differentially expressed in the two comparison groups, then constructed the competing endogenous RNA (ceRNA) network, which showed up-regulated miRNAs with decreased expression of protein-coding and circRNA transcripts, while the down-regulated miRNAs showed the reverse results. There were 190 nodes (Fig 5), which suggested that the ceRNA network was very dense.
Subsequently, a new hypoxic ceRNA network were constructed. It is worth mentioning that the 2 lncRNAs (LOC100335190 and TCONS_00053712) were the differentially expressed lncRNAs in LAC-vs-LWQY comparison group. In the hypoxic ceRNA network, 9 novel circRNA transcripts were predicted to interact with MAPKAPK3, NFATC2, ATP7A and DIAPH1 through miR-195-x, which indicated that miR-195-x may play an important role in hypoxic adaptation for yak in plateau.

Discussion
High-altitude adaptation is an extremely complex biological process that involves millions of molecules and biological pathways. A large number of SNPs and copy number variations (CNVs) have been widely researched and showed that many gene maybe the targets for further high-altitude adaptation mechanism studies in pig [40] and human [41], and accumulating evidence has associated small non-coding RNAs, long non-coding RNAs and circular RNAs with hypoxia signaling pathway in different cell types. Yak provide an ideal model for deciphering the mechanism of high-altitude adaptation for their well adaptation in the harsh and extreme environment. Despite the identification of SNPs, CNVs [42], mRNA and miRNA [13] expression profiles for adaptation of yak, the whole transcriptome and the connection between mRNA and non-coding RNAs remains poorly understood.
Comparing the whole transcriptomes of yak and cattle can provide insight into the regulatory mechanisms of high-altitude adaptation at a comprehensive view, including among mRNAs, miRNAs, lncRNAs, and circRNAs. Here, we performed whole transcriptomes analysis in the heart of Leiwuqi (LWQ) yak and their relatives, the migration cattle in LWQ yak habitat (HAC) and low-altitude cattle (LAC), respectively, and identified comparable numbers of transcriptomes species using yak and cattle genome, respectively. In addition, we also constructed the hypoxia-adaptation ceRNA network. PIK3R5 is one subunit of phosphoinositide 3-kinaseγ, which is a lipid kinase that belongs to class IB subdivision of Phosphoinositide 3-kinases (PI3Ks). Previous studies show that the activated PI3K-Akt pathway in hypoxic endothelial cells is essential for cell survival and angiogenesis, and involved in hypoxia-inducible factor 1αsignaling pathway [43][44]. Hence, we speculate that the activated PIK3R5 in heart of yak may play a key role in protecting cells from hypoxia via PI3K-AKT signaling pathway.
NFAT transcription factors play critical roles in gene transcription during immune responses, as two most prominent NFAT family members, NFATC1 and NFATC2 plays an important role in controlling both T and B cell activation and differentiation [45][46]. In addition, our results also showed that DEGs were not only associated with above mentioned pathways but also with other pathways, for instance, the VEGF signaling pathway, and cardiac muscle contraction. Taken together, these results indicated the potential roles of these DEGs in high altitude adaptation through activating important pathways. lncRNA is one major class of ncRNAs and could act as important regulators of gene expression in a wide variety of biological processes, including energy balance and immune responses. It has been reported that the energy stress-induced lincRNA-FILNC1 could represses c-Myc-mediated energy metabolism and lincRNA-Cox2 could be induced in mouse macrophages after activation of Toll-like receptors so as to detect microbes and alter the immune system to respond [47]. This experiment is the first to show expression profiles of lncRNA transcripts in yak heart and differential expressed lncRNA transcripts in response to high-altitude using whole transcriptome sequencing. Through functional analysis, 168 target genes of differential expressed lncRNA transcripts were enriched in 388 GO terms and 17 signaling pathways, the results of which demonstrate that lncRNAs have important roles in high-altitude animals. Differ from protein-coding transcripts, although we found substantial differentially expressed lncRNAs that may be associated with high altitude adaptation, no co-differentially expressed lncRNAs were identified between LAC-vs-LWQY and LAC-vs-HAC groups.
Such discrepancy results could be partly explained by the lncRNAs are often less conserved at the level of primary sequence and hinders the discovery of lncRNA orthologs in yak and cattle genomes by sequence homology, and since diverging five million years ago, yak and cattle underwent distinct selective pressure resulting in inherent genetic difference which may to some extent disturb the results. Although data from whole transcriptome of yak and its relatives didn't reveal the codifferential expressed lncRNAs, but the link between differential expression lncRNAs in LAC-vs-LWQY and LAC-vs-HAC groups and hypoxia adaptation appears to be biologically relevant.
In this study, we identified a total of 58 co-differentially expressed miRNAs between LAC-vs-LWQY and LAC-vs-HAC groups, more than that reported by a previous study (29 miRNAs) [13], which compared the miRNAs expression profiles in heart tissues between yak and cattle. The previous study showed that target genes of miRNAs from heart were mainly associated with hypoxia-related functions [13], including the HIFa signaling pathway, PI3K-Akt signaling pathway, and p53 signaling pathway. However, we found that the co-differentially target genes of the co-differentially expressed miRNAs were significantly involved in the VEGF signaling pathway, cAMP signaling pathway, T cell receptor signaling pathway, and B cell receptor signaling pathway, which are involved in energy metabolism, immunity, and angiogenesis. Interestingly, we found that the targets were also involved in longevity regulating pathway, osteoclast differentiation and etc. Hence, we suspected that the codifferentially expressed miRNAs may be involved in the regulation of not only immune-, energy-, and angiogenesis responses but also other biology processes.
In addition to mRNAs, lncRNAs, and miRNAs, circRNAs are a distinct class of newly discovered endogenous ncRNAs. Indeed, circRNAs were found to play crucial roles in various biological and pathological processes across a wide range of organisms [48][49]. However, there is little known about the features of circRNAs in yak. To explore the expression profiles of circRNAs and their potential function in the regulation in high-altitude adaptation, we first compared the differentially expressed circRNAs between yak and its relatives. Finally, we found that a total of 53 circRNAs showed codifferential expression between LAC-vs-LWQY and LAC-vs-HAC comparison groups. An analysis of circRNA-miRNA associations showed that 53 co-differentially expressed circRNAs can regulate 57 identified co-differentially expressed miRNAs, suggesting circRNAs may also play vital roles in the regulation of a wide range of signaling pathway hence involved in the regulation mechanism of highaltitude adaptation.
Integration of multi-omics can generate new information and improve classification accuracy that are not accessible by analysis of single datasets alone [50][51]. For instance, Liu et al. demonstrated that the integration of multiple omics improved the validity of functional analysis of differentially expressed miRNAs using intersection genes and constructed a ceRNA relationship network in broody chickens [52]. Drake et al. integrated genomic, transcriptomic, and phosphoproteomic datasets to reveal a map of activated signaling pathways in castration-resistant prostate cancer [53]. In the present study, we integrated multiple omics to identify altered genes and processes in high-altitude adaptation in Bovidae, and constructed a ceRNA relationship network to provide a reliable reference for further investigation.

Conclusion
Taken together, this study has investigated whole transcriptome profiles of nine bovidae hearts from yak and their relatives (the migration cattle in yak habitat and low-altitude cattle). We have identified 178 protein-coding transcripts, 58 miRNAs, and 53 circRNAs were differentially co-expressed in LACvs-LWQY and LAC-vs-HAC comparison groups, and In silico analysis indicated that those transcripts can regulate a wide range of biological processes in high-altitude adaptation in Bovidae. More interestingly, the present study indicated that 8 miRNAs and 15 circRNAs were directly or indirectly interact with 6 protein-coding transcripts (MAPKAPK3, PXN, NFATC2, ATP7A, DIAPH1 and F2R) in the ceRNA network. Our data reinforce the view that the molecular network of high-altitude adaptation is composed of several protein-coding transcripts and non-coding transcripts that involved in multiply regulatory interactions. The detailed mechanism for how this network crosstalk each other in heart of bovidae are worthy of future research efforts.

Tissue Samples and RNA Extraction
A total of 9 healthy unrelated adult females for Leiwuqi (LWQ) yak and their relatives, the migration cattle in LWQ yak habitat (altitude, 3600m, HAC) and low-altitude cattle (altitude, 480m, LAC), were randomly selected for heart tissue sampling, each group includes three healthy individuals of similar body weight. The LWQ yak and the HAC cattle were reared under similar environmental and feeding conditions in Leiwuqi country (Changdu, Tibet, China), while the LAC cattle were reared under the same feeding conditions as LWQ yak and HAC cattle in Wenchuan country (Aba prefecture, Sichuan, China). With the Trizol reagent (Invitrogen, USA), total RNA was extracted according to the manufacturer's instructions. The RNA concentration, purity and integrity were determined by NanoDrop 2000 spectrophotometer (Thermo Scientific).
Next generation sequencing and data processing The next generation sequencing (NGS) service was provided by Gene Denovo Biotechnology Co.
(Guangzhou, China). For miRNA sequencing, the miRNA library was prepared according to experimental procedures reported previously [13], and was sequenced using Illumina HiSeq TM 2500 following the vendor's recommended protocol. Raw data were successively filtered according to the following rules:1) Removing low quality reads containing more than one low quality (Q-value≤20) base or containing unknown nucleotides (N); 2) removing reads containing 5' adaptor, no 3' adaptor or no insertion sequence; 3) removing reads containing poly (A) in small RNA fragment or shorter than 18 nt (not include adaptor); 4) Bovine known classes of RNAs (mRNA, rRNA, tRNA, snRNA, scRNA, snoRNA and repeats) were subsequently removed through searching against three databases, including GenBank database (Realse 209.0), Rfam and reference Genome. For cattle, All the retains reads were then searched against miRBase database (Release 22) to identify exist miRNAs and known miRNAs, then the unannotated reads were aligned with cattle genome to identify novel miRNAs. For yak, the retains reads were mapped to the known miRNAs in miRBase 22.0 to obtain known miRNAs, and the unannotated reads were aligned with yak genome to obtain novel miRNAs.
The cDNA libraries of RNA (mRNA/lncRNA/circRNA) were generated using the mRNA-Seq sample preparation kit (Illumina) after ribosomal RNA was removed by the Epicentre Ribo-Zero rRNA kit (Epicentre, USA), and was sequenced using Illumina HiSeq TM 4000 following the vendor's recommended protocol. The detail data processing protocols are as reported previously [54]. In brief, the reads with adaptor contamination, low quality bases and undetermined bases were removed using the Cutadapt. Bowtie2 and Tophat2 were used to map the reads to the bos grunniens or bos taurus genome. Putative protein-coding RNAs were filtered out using a minimum length and exon number threshold. Transcripts of more than 200 nt and more than two exons were selected as lncRNA candidates and further analysis. The programs Coding-Non-coding-Index (version 2) [55], Coding Potential Calculator (CPC, http://cpc2.cbi.pku.edu.cn/) [56] and phylogenetic codon substitution frequency (PhyloCSF) [57] were used to predict the protein-coding potential of new transcripts with default parameters. The intersection of the results without protein-coding potential was yielded lncRNA transcripts.
Firstly, CIRI employed BWA to align the clean reads with a reference genome to generate SAM files.
Then, the CIGAR values and the junction reads with paired chiastic clipping (PCC) signals were analyzed and detected by using the CIRI software. Secondly, the SAM files were scanned by CIRI to detect additional junction reads exclude false positive candidates. Differentially expressed circRNAs were identified using EBSeq.
The expression of all the transcripts was estimated with the Stringtie and Ballgown software [59].
Furthermore, the edger package (http://www.bioconductor.org/packages/release/bioc/html/edgeR.html) was used to identify differentially expression transcripts across samples or groups, we identified mRNA and lncRNA with a false discovery rate (FDR) 0.05 and |log2(Fold Change) |≥1 in a comparison as significant differentially expression transcripts (DETs), and miRNA and circRNA with |log2(Fold Change) |≥1 and P 0.05. It's worth noting that one cattle gene ID corresponds to one or more yak gene ID, and the codifferentially transcripts were analyzed based on cattle gene ID.

The association analysis between lncRNA and mRNA
Previous studies point towards an increasing number of lncRNAs function in specific physiological and pathological contexts are associated with mRNA expression level [23], to reveal the association between the lncRNA and mRNA, the antisense-lncRNA, cis-lncRNA and trans-lncRNA were predicted as follows: (1)  MiRNA target prediction and Construction of ceRNA network Three softwares, mireap, miRanda and TargetScan, were used to predict the miRNAs and circRNAs targets. MiRNA sequences and family information were obtained from TargetScan website (http://www.targetscan.org/).
The ceRNA network was constructed based on ceRNA theory [26,33] as follows: (1) Expression correlation between mRNA-miRNA or lncRNA-miRNA was evaluated using the Spearman Rank correlation coefficient (SCC), Pairs with SCC -0.7 were selected as co-expressed negatively lncRNA-miRNA pairs or mRNA-miRNA pairs. Both mRNA and lncRNA were miRNA targets genes, and all RNAs were differentially expressed. (2) Expressed correlation between lncRNA-mRNA was evaluated using the Pearson correlation coefficient, pairs with Pearson correlation coefficient 0.9 were selected as coexpressed lncRNA-mRNA pairs, both mRNA and lncRNA in this pair were targeted and co-expressed negatively with a common miRNA. (3) Hypergeometric cumulative distribution function test was used to test whether the common miRNA sponges between the two genes were significant. As a result, only the gene pairs with a p-value less than 0.05 were selected. The lncRNA-miRNA-mRNA network was constructed by assembling all co-expression competing triplets, which were identified above, and was visualized using Cytoscape software (v3.6.0, http://www.cytoscape.org/).     The KEGG pathway analysis of the 125 differentially co-expressed genes targeted by the 58 differentially co-expressed miRNAs. The predicted hypoxic ceRNA network. The red square, green circle, yellow triangle and dark blue rectangle indicates protein-coding transcript, circRNA transcript, miRNA and lncRNA transcript, respectively.

Supplementary Files
This is a list of supplementary files associated with this preprint. Click to download.