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

doi:10.21203/rs.2.10149/v1

Download PDF

Research article

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

https://doi.org/10.21203/rs.2.10149/v1

This work is licensed under a CC BY 4.0 License

Version 1

posted

You are reading this latest preprint version

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-HAC，respectively, but no lncRNA transcripts co-differentially expressed in those two comparison groups. Among the 58 miRNAs that were co-differentially expressed, 18 were up-regulated 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, ATP7A, DIAPH1 and F2R genes. Conclusion: In conclusion, 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 hypoxic adaptation regulatory pathways. The detailed mechanism for how this network crosstalk each other in heart of Bovidae are worthy of future research efforts.

Epigenetics & Genomics

Yak

Cattle

Whole-transcriptome

Co-differentially expressed transcripts

Hypoxic adaptation

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 hypoxia-inducible 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 (Supplementary file1). Paired-end sequencing of the miRNA library yielded 145,480,110 raw reads and 132,141,703 clean reads from the nine libraries, with an average of 14,682,411 clean reads in each library (Supplementary file2).

The expression analysis and function annotation of protein-coding transcripts

In the strand-specific sequencing libraries, we identified 63,234 protein-coding transcripts including 44,233 known and 7,957 novel protein-coding transcripts in cattle, while 18,891 known and 8,381 novel protein-coding transcripts in yak (Supplementary file3, Table 1). After transcripts abundances were quantified by Stringtie and Ballgown software, the average expression level of novel protein-coding transcripts (2.93) is about twenty-eight percent of known protein-coding transcripts (10.37) in cattle, and the average expression level of novel protein-coding transcripts (7.43) is about thirty-two percent of known protein-coding transcripts (23.48) in LWQY (Supplementary file3). 27,222, 48,767 and 33,235 protein-coding transcripts were detected in the whole transcripts in LWQY, HAC, and LAC groups, respectively (Table 1), and 32,094 protein-coding transcripts of yak were mapped to 23,346 cattle transcripts (Supplementary file3). A total of 18,437 transcripts of yak and 14,297 protein-coding transcripts of cattle were co-expressed in LWQY-vs-HAC-vs-LAC comparison groups. A total of 2,351, 1,452 and 185 protein-coding transcripts were specifically expressed in LWQY, HAC and LAC, respectively (Supplementary file3). Furthermore, a large number of alternative splicing events were 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.

Overall, 5,637 (including 3,215 decreased and 2,422 increased) and 230 (including 191 decreased and 39 increased) significantly DETs were detected in the LAC-vs-LWQY (Supplementary file4) and LAC-vs-HAC group (Supplementary file5), respectively. While, 178 significantly DETs including 6 down-regulated and 172 up-regulated DETs were co-differentially discovered in both LAC-vs-LWQY and LAC-vs-HAC groups (Fig1B). GO enrichment (Supplementary file6) and KEGG pathway analysis (Supplementary file7 and 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 with immune system and protein metabolism in mammals.

The expression analysis and function annotation of lncRNA transcripts

In total, there were 4,657 lncRNA transcripts. 1429, 3556 and 1540 lncRNA transcripts were detected in the whole transcriptome of LWQY, HAC, and LAC, respectively (Supplementary file8). The lncRNA transcripts of LWQY corresponded to 1222 lncRNA genes, while the lncRNA transcripts of HAC and LAC corresponded to 2874 lncRNA genes, and there were 864 lncRNA transcripts (691 lncRNA transcripts genes) of yak mapped to 689 lncRNA transcripts of cattle (620 lncRNA genes). Furthermore, in LWQY-vs-HAC-vs-LAC comparison group, 420 lncRNA transcripts of yak and 323 lncRNA transcripts of cattle were co-expressed, respectively. This difference was caused by that one lncRNA transcripts of cattle corresponds to one or more lncRNA transcripts of yak. In addition, 131, 11 and 3 lncRNA transcripts were specific expressed in LWQY, HAC, and LAC group, respectively (Supplementary file9). The specific expressed lncRNAs in different groups may be a reflection of its physiological role in the regulation of hypoxic adaptation in Bovidae.

A boxplot of average expression level, histogram of lengths, proportion of exon number per transcript, and GC content for protein-coding and lncRNA transcripts are shown in Fig 3. 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 protein-coding 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-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 protein-coding genes of the differential lncRNA transcripts. Functional analysis showed that these target genes were significantly enriched in 388 GO terms (304 under biological process, 32 under cellular component, and 52 under molecular function), and many terms were related to morphogenesis, immune response and gene expression (Supplementary Fig1). For example, the biological process involved cell morphogenesis and cell morphogenesis involved in differentiation, immune system development, regulation of cell differentiation, methylation and histone deacetylation. In addition, the target genes were enriched in 17 pathways, several of which were associated with environmental adaptation, immune system, and carbohydrate metabolism pathways which are very important for yak hypoxia adaptation (Supplementary Fig2). These results indicated that lncRNAs took part in several biological processes in the yak heart during hypoxia adaptation.

The expression analysis and function annotation of miRNA transcripts

A total of 1937 miRNAs were identified in the nine small RNA libraries (Supplementary file10), which length were distributed in 18-24 nucleotides. Total 831, 1713 and 1713 miRNAs were detected from LWQY, HAC and LAC, respectively (Supplementary file10). The top 10 mature miRNAs with the highest expression comprised more than 61% of all miRNAs in LWQY, HAC and LAC, showing a relatively abundant distribution. Furthermore, miR-99-x were the general abundant miRNA overall, together accounting for 20%, 8%, and 4% of the total in LWQY, HAC and LAC, respectively (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-miRNA-pathway 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 (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.

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 (Supplementary file17) and LAC-vs-HAC (Supplementary file18) comparison groups, respectively, and 53 up-regulated differential circRNAs 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.

Based on the combining analysis of the ceRNA network, 8 miRNAs and 15 circRNAs were predicted to directly or indirectly interact with 6 mRNAs (MAPKAPK3 [36], PXN [37], NFATC2 [38], ATP7A [39], DIAPH1 and F2R), which were reported to involved in hypoxia adaptation related signaling pathways. 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.

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.

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 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 co-differential 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 co-differentially 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 co-differential 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 high-altitude 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.

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 LAC-vs-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.

CIRI software (https://sourceforge.net/projects/ciri/) was used to identify circular RNAs (circRNAs) [58]. 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 co-differentially 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) LncRNA is involved in the regulation of many post-transcriptional processes, similar to small RNAs such as miRNA and snoRNA, which are often associated with complementary pairing of bases. A portion of antisense lncRNA may regulate gene silencing, transcription, and mRNA stability via binding to the sense strand of mRNA. To reveal the interaction between antisense lncRNA and mRNA, RNAplex software was used to predict the complementary bases between antisense lncRNA and mRNA. (2) The lncRNAs annotated in the "unknown region" were used for cis regulation analysis. The lncRNAs located within 10 kb upstream or downstream of a gene may overlap with the region of the cis-acting element and participate in mRNA transcriptional regulation, we named these lncRNAs as cis-lncRNAs. (3) The principle of trans target gene prediction of lncRNA is that the function of lncRNA is not related with the location of the coding gene, but co-expression with the protein coding gene. The co-expression analysis of lncRNA and protein-coding gene was used for predicting the target genes between samples. LncTar was used to calculate the free energy between lnRNAs and mRNAs and to predict the targets of lncRNAs [60].

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 co-expressed 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.

P-value=1-F(x/U,M,N)

For a given each gene pair (A, B), we denoted all their regulator miRNAs as miRNA set C (regulating gene A) and D (regulating gene B). Where X stands for the number of common miRNAs that regulate both of the two genes, U means the number of all the miRNAs in this work, M presents the size of miRNA set C and N present the size of miRNA set D.

Functional enrichment and Visualization of ceRNA network

To annotate the function of DEGs, Gene Ontology (GO) enrichment and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analyses were performed using Cytoscape.

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

LWQ: Leiwuqi yak; HAC: the migration catlle in LWQ yak habitat; LAC: low-altitude cattle; CETs: co-differentially expressed protein-coding transcripts; ceRNA: competing endogenous RNA; ncRNA: non-coding RNA; lncRNA: long noncoding RNA; circRNA: circular RNA; miRNA: microRNA; 3’-UTR: 3’-untranslated region; HIF1α: hypoxia-inducible factors 1α; MREs: miRNA response elements; NGS: next generation sequencing; GO: Gene Ontology; KEGG: Kyoto Encyclopedia of Genes and Genomes; COX7C: cytochrome c oxidase subunit 7C; MAPKAPK3: mitogen-activated protein kinase-activated protein kinase 3; BAD: BCL2 associated agonist of cell death; VAV1: vav guanine nucleotide exchange factor 1; NFATC2: nuclear factor of activated T cells 2; NFATC1: nuclear factor of activated T cells 1; PIK3R5: phosphoinositide-3-kinase regulatory subunit 5; PXN: paxillin; ATP7A: ATPase copper transporting alpha; DIAPH1: diaphanous related formin 1; F2R: coagulation factor II thrombin receptor; CNV: copy number variations.

Acknowledgements

Not applicable

Funding

This study was supported by the Key Laboratory of Qinghai-Tibetan Plateau Animal Genetic Resource Reservation and Utilization Open Foundation in designing the study and sample collection; the Open Project Program of State Key Laboratory of Barley and Yak Germplasm Resources and Genetic Improvement, Tibet Academy of Agricultural and Animal Husbandry Sciences (XZNKY-2019-C-007K08) in analysis and interpretation of the data; and the Program of National Beef Cattle and Yak Industrial Technology System (CARS-37) in writing the manuscript.

Availability of data and materials

The datasets supporting the conclusions of this article are included as supplementary tables: Supplementary file1-19. The raw data would be available by contact with the corresponding author upon reasonable request.

Authors’ contributions

WH, ZJC, and CZX conceived and designed the experiments, WH, CZ and ZJJ analyzed the data, WJK, ZCF, CHW and XJW collected the samples. WH wrote the paper. All authors read and approved the final manuscript.

Ethics approval and consent to participate

All animal experiments were performed according to the regulations for the Administration of Affairs Concerning Experimental Animals (Ministry of Science and Technology, China, revised in June 2004) and approved by the Institution Animal Care and Use Committee in the Southwest Minzu University, Chengdu, China. We obtained written permission from the owner of the yak and cattle to take samples.

Consent for Publication

Not applicable.

Competing interests

The authors declared that they have no competing interests.

1. Beja-Pereira A, Caramelli D, Lalueza-Fox C, Vernesi C, Ferrand N, Casoli A, et al. The origin of European cattle: evidence from modern and ancient DNA[J]. Proceedings of the National Academy of Sciences. 2006; 103(21): 8113-8118.

2. Qiu Q, Zhang G, Ma T, Qian W, Wang J, Ye Z, et al. The yak genome and adaptation to life at high altitude[J]. Nature genetics. 2012; 44(8): 946-949.

3. Wiener G, Han J L, Long R J. The yak. Second edition. Bangkok: FAO Regional Office for Asia and the Pacific Food and Agriculture Organization of the United Nations. 2003

4. Weir EK, Tucker A, Reeves JT, Will DH, Grover RF. The genetic factor influencing pulmonary hypertension in cattle at high altitude[J]. Cardiovascular Research. 1974; 8(6): 745-749.

5. Hecht HH, Kuida H, Lange RL, Thorme JL, Browen AM. Brisket disease: II. Clinical features and hemodynamic observations in altitude-dependent right heart failure of cattle[J]. The American journal of medicine. 1962; 32 (2): 171-183.

6. Shao B, Long R, Ding Y, Wang J, Ding L, Wang H. Morphological adaptations of yak (Bos grunniens) tongue to the foraging environment of the Qinghai-Tibetan Plateau[J]. Journal of Animal Science. 2010; 88(8): 2594-2603.

7. Wang H, Long R, Liang JB, Guo XS, Ding LM, Shang ZH. Comparison of nitrogen metabolism in Yak (Bos grunniens) and indigenous cattle (Bos taurus) on the Qinghai-Tibetan Plateau[J]. Asian-Australasian Journal of Animal Sciences. 2011; 24(6): 766-773.

8. Hoppeler H, Kleinert E, Schlegel C, Claassen H, Howald H, Kayar1 SR, et al. II. Morphological adaptations of human skeletal muscle to chronic hypoxia[J]. International journal of sports medicine. 1990; 11(S1): S3-S9.

9. Ahmad K S, Hameed M, Fatima S, Ashraf M, Ahmad F, Naseer M, et al. Morpho-anatomical and physiological adaptations to high altitude in some Aveneae grasses from Neelum Valley[J]. Western Himalayan Kashmir, Acta physiologiae plantarum. 2016; 38(4): 93.

10. Monge C, Leon-Velarde F. Physiological adaptation to high altitude: oxygen transport in mammals and birds[J]. Physiological Reviews. 1991; 71(4): 1135-1172.

11. Wang K, Yang Y, Wang L, Ma T, Shang H, Ding L, et al. Different gene expressions between cattle and yak provide insights into high‐altitude adaptation[J]. Animal genetics. 2016; 47(1): 28-35.

12. Qi X, Zhang Q, He Y, Yang L, Zhang X, Shi P, et al. The transcriptomic landscape of yaks reveals molecular pathways for high altitude adaptation[J]. Genome biology and evolution. 2019; 11(1): 72-85.

13. Guan J, Long K, Ma J, Zhang J, He D, Jin L, et al. Comparative analysis of the microRNA transcriptome between yak and cattle provides insight into high-altitude adaptation[J]. PeerJ. 2017; 5: e3959.

14. Liu L, Xiao Q, Gilbert ER, Cui Z, Zhao, Wang Y, et al. Whole-transcriptome analysis of atrophic ovaries in broody chickens reveals regulatory pathways associated with proliferation and apoptosis[J]. Scientific reports. 2018; 8(1): 7231.

15. He JH, Han ZP, Zou MX, Wang L, Lv YB, Zhou JB, et al. Analyzing the LncRNA, miRNA, and mRNA regulatory network in prostate cancer with bioinformatics software[J]. Journal of Computational Biology. 2018; 25(2): 146-157.

16. Chen B, Huang S. Circular RNA: an emerging non-coding RNA as a regulator and biomarker in cancer[J]. Cancer letters. 2018; 418: 41-50.

17. Shih JW, Chiang WF, Wu A, Wu MH, Wang LY, Yu YL, et al. Long noncoding RNA LncHIFCAR/MIR31HG is a HIF-1α co-activator driving oral cancer progression[J]. Nature communications. 2017; 8: 15874.

18. Choudhry H, Harris A L. Advances in hypoxia-inducible factor biology[J]. Cell metabolism. 2018; 27(2): 281-298.

19. Massart J, Sjögren RJO, Lundell LS, Leonidas SL, Jonathan MM, Niclas F, et al. Altered miRNA-29 Expression in Type 2 Diabetes Influences Glucose and Lipid Metabolism in Skeletal Muscle[J]. Diabetes. 2017; 66(7):1807-1818.

20. Goyal N, Kesharwani D, Datta M. Lnc-ing non-coding RNAs with metabolism and diabetes: roles of lncRNAs[J]. Cellular and Molecular Life Sciences. 2018; 75(10): 1827-1837.

21. Anastasiadou E, Jacob L S, Slack F J. Non-coding RNA networks in cancer[J]. Nature Reviews Cancer. 2018; 18(1): 5-18.

22. Filipowicz W, Bhattacharyya SN, Sonenberg N. Mechanisms of post-transcriptional regulation by microRNAs: are the answers in sight? [J]. Nature reviews genetics. 2008; 9(2): 102-114.

23. Wang KC, Chang HY. Molecular mechanisms of long noncoding RNAs[J]. Molecular cell. 2011; 43(6): 904-914.

24. Bhan A, Deb P, Shihabeddin N, Ansari KI, Brotto M, Mandal SS. Histone methylase MLL1 coordinates with HIF and regulate lncRNA HOTAIR expression under hypoxia[J]. Gene. 2017; 629: 16-28.

25. Hansen TB, Jensen TI, Clausen BH, Bramsen JB, Finsen B, Damgaard CK, et al. Natural RNA circles function as efficient microRNA sponges[J]. Nature. 2013; 495(7441): 384-388. Doi: 10.1038/nature11993

26. Salmena L, Poliseno L, Tay Y, Pandolfi PP. A ceRNA hypothesis: the Rosetta Stone of a hidden RNA language? [J]. Cell. 2011; 146(3): 353-358.

27. Tong Q, Zheng L, Lin L, Li B, Wang D, Huang C, et al. VEGF is upregulated by hypoxia-induced mitogenic factor via the PI-3K/Akt-NF-κB signaling pathway[J]. Respiratory research. 2006; 7(1): 37.

28. Knerr I, Schubert SW, Wich C, Amann K, Aigner T, Vogler T, et al. Stimulation of GCMa and syncytin via cAMP mediated PKA signaling in human trophoblastic cells under normoxic and hypoxic conditions[J]. FEBS letters. 2005; 579(18): 3991-3998.

29. Vaupel P, Multhoff G. Hypoxia-/HIF-1α-Driven Factors of the Tumor Microenvironment Impeding Antitumor Immune Responses and Promoting Malignant Progression//Oxygen Transport to Tissue XL[J]. Springer, Cham. 2018; 171-175.

30. Amoasii L, Olson E N, Bassel-Duby R. Control of muscle metabolism by the mediator complex[J]. Cold Spring Harbor perspectives in medicine. 2018; 8(2): a029843.

31. Klecha AJ, Genaro AM, Gorelik G, Barreiro Arcos ML, Silberman DM, Schuman M, et al. Integrative study of hypothalamus–pituitary–thyroid–immune system interaction: thyroid hormone-mediated modulation of lymphocyte activity through the protein kinase C signaling pathway[J]. Journal of Endocrinology. 2006; 189(1): 45-55.

32. Häussinger D, Schliess F. Glutamine metabolism and signaling in the liver[J]. Front Biosci. 2007; 12(4): 371-391.

33. Xu X, Zhou X, Wang R, Peng W, An Y, Chen L. Functional analysis of long intergenic non-coding RNAs in phosphate-starved rice using competing endogenous RNA network[J]. Scientific reports. 2016; 6: 20715.

34. Lan X, Wang Y, Tian K, Ye F, Yin H, Zhao X, et al. Integrated host and viral transcriptome analyses reveal pathology and inflammatory response mechanisms to ALV-J injection in SPF chickens[J]. Scientific Reports. 2017; 7: 46156.

35. Han D, Zhang Y, Chen J, Hua G, Li J, Deng X, Deng X. Transcriptome analyses of differential gene expression in the bursa of Fabricius between Silky Fowl and White Leghorn[J]. Scientific reports. 2017; 7: 45959.

36. Mann S, Sipka A S, Grenier J K. The degree of postpartum metabolic challenge in dairy cows is associated with peripheral blood mononuclear cell transcriptome changes of the innate immune system[J]. Developmental & Comparative Immunology. 2019; 93: 28-36.

37. Veith C, Schermuly RT, Brandes RP, Weissmann N. Molecular mechanisms of hypoxia‐inducible factor‐induced pulmonary arterial smooth muscle cell alterations in pulmonary hypertension[J]. The Journal of physiology. 2016; 594(5): 1167-1177.

38. Peng S L, Gerth A J, Ranger A M, Glimcher LH. NFATc1 and NFATc2 together control both T and B cell activation and differentiation[J]. Immunity. 2001; 14(1): 13-20.

39. Xie L, Collins J F. Transcriptional regulation of the Menkes copper ATPase (Atp7a) gene by hypoxia-inducible factor (HIF2α) in intestinal epithelial cells[J]. American Journal of Physiology-Cell Physiology. 2011; 300(6): C1298-C1305.

40. Tang Q, Gu Y, Zhou X, Jin L, Guan J, Liu R, Li J, Long K, Tian S, Che T, Hu S, Liang Y, Yang X, Tao X, Zhong Z, Wang G, Chen X, Li D, Ma J, Wang X, Mai M, Jiang A, Luo X, LV X, Gladyshev V, Li X, Li M. Comparative transcriptomics of 5 high-altitude vertebrates and their low-altitude relatives[J]. GigaScience. 2017; 6(12): 1-9.

41. Yang J, Jin Z, Chen J, Huang X, Li X, Liang Y, Mao J, Chen X, Zheng Z, Bakshi A, Zheng D, Zheng M, Wray N, Visscher P, Lu F, Qu J. Genetic signatures of high-altitude adaptation in Tibetans[J]. Proceedings of the National Academy of Sciences. 2017; 114(16): 4189-4194.

42. Wang H, Chai Z, Hu D, Ji Q, Xin J, Zhang C, Zhong J. A global analysis of CNVs in diverse yak populations using whole-genome resequencing[J]. BMC genomics. 2019; 20(1): 61.

43. Hung S, Pochampally R, Chen S, Hsu S, Prockop D. Angiogenic effects of human multipotent stromal cell conditioned medium activate the PI3K‐Akt pathway in hypoxic endothelial cells to inhibit apoptosis, increase survival, and stimulate angiogenesis[J]. Stem cells, 2007, 25(9): 2363-2370.

44. Zhou J, Schmid T, Frank R, Brüne B. PI3K/Akt is required for heat shock proteins to protect hypoxia-inducible factor 1α from pVHL-independent degradation[J]. Journal of Biological Chemistry, 2004, 279(14): 13506-13513. 52. Macian F. NFAT proteins: key regulators of T-cell development and function[J]. Nature Reviews Immunology, 2005, 5(6): 472-484.

45. Macian F. NFAT proteins: key regulators of T-cell development and function[J]. Nature Reviews Immunology, 2005, 5(6): 472-484.

46. Giampaolo S, Wójcik G, Klein-Hessling S, Serfling E, Patra A. B cell development is critically dependent on NFATc1 activity[J]. Cellular & molecular immunology, 2018: 1.

47. Xiao Z, Han L, Lee H, Zhuang L, Zhang Y, Baddour J, Nagrath D, Wood C, Gu J, Wu X, Liang H, Gan B. Energy stress-induced lncRNA FILNC1 represses c-Myc-mediated energy metabolism and inhibits renal tumor development[J]. Nature Communications, 2017, 8(1):783.

48. Dai X, Chen C, Yang Q,Xue J, Chen X, Sun B, Luo F, Liu X, Xiao T, Xu H, Sun Q, Zhang A, Liu Q. Exosomal circRNA_100284 from arsenite-transformed cells, via microRNA-217 regulation of EZH2, is involved in the malignant transformation of human hepatic cells by accelerating the cell cycle and promoting cell proliferation[J]. Cell Death & Disease, 2018, 9(5):454.

49. Xu S, Zhou L, Ponnusamy M, Zhang L, Dong Y, Zhang Y, Wang Q, Liu J, Wang K. A comprehensive review of circRNA: from purification and identification to disease marker potential [J]. Peer J, 2018, 6: e5503.

50. Zhang W, Li F, Nie L. Integrating multiple 'omics' analysis for microbial biology: application and methodologies[J]. Microbiology, 2010, 156(Pt 2):287-301

51. Vilne B, Schunkert H. Integrating Genes Affecting Coronary Artery Disease in Functional Networks by Multi-OMICs Approach[J]. Frontiers in Cardiovascular Medicine, 2018, 5: 89.

52. Liu L, Xiao Q, Gilbert E, Cui Z, Zhao X, Wang Y, Yin H, Li D, Zhang H, Zhu Q. Whole-transcriptome analysis of atrophic ovaries in broody chickens reveals regulatory pathways associated with proliferation and apoptosis[J]. Sci Rep, 2018, 8(1): 7231.

53. Drake J, Paull E, Graham N, Lee J, Smith B, Titz B, Stoyanova T, Faltermeier C, Uzunangelov V, Carlin D, Fleming D, Wong C, Newton Y, Sudha S, Vashisht A, Huang J, Wohlschlegel J, Graeber T, Witte O, Stuart J. Phosphoproteome Integration Reveals Patient-Specific Networks in Prostate Cancer[J]. Cell, 2016, 166: 1041-1054.

54. Ye P, Shi Y, An N, Zhou Q, Guo J, Long X. miR-145 overexpression triggers alteration of the whole transcriptome and inhibits breast cancer development[J]. Biomedicine & Pharmacotherapy. 2018; 100: 72-82.

55. Sun L, Luo H, Bu D, Zhao G, Yu K, Zhang C, et al. Utilizing sequence intrinsic composition to classify protein-coding and long non-coding transcripts[J]. Nucleic acids research. 2013; 41(17): e166.

56. Kong L, Zhang Y, Ye ZQ, Liu XQ, Zhao SQ, Wei L, et al. CPC: assess the protein-coding potential of transcripts using sequence features and support vector machine[J]. Nucleic acids research. 2007; 35(suppl_2): W345-W349.

57. Lin M F, Jungreis I, Kellis M. PhyloCSF: a comparative genomics method to distinguish protein coding and non-coding regions[J]. Bioinformatics. 2011; 27(13): i275-i282.

58. Gao Y, Wang J, Zhao F. CIRI: an efficient and unbiased algorithm for de novo circular RNA identification[J]. Genome biology. 2015; 16(1): 4.

59. Pertea M, Kim D, Pertea G M, Leek JT, Salzberg SL. Transcript-level expression analysis of RNA-seq experiments with HISAT, StringTie and Ballgown[J]. Nature protocols. 2016; 11(9): 1650-1667.

60. Li J, Ma W, Zeng P, Wang J, Geng B, Yang J. LncTar: a tool for predicting the RNA targets of long noncoding RNAs[J]. Briefings in bioinformatics. 2014; 16(5): 806-812.

Due to technical limitations, all Tables are only available as a download in the supplemental files section.

Download PDF

Version 1

posted

You are reading this latest preprint version

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

Status:

Version 1

Abstract

Figures

Background

Results

Overview of RNA-sequencing

The expression analysis and function annotation of protein-coding transcripts

The expression analysis and function annotation of lncRNA transcripts

The expression analysis and function annotation of miRNA transcripts

The circRNAs expression analysis and circRNA-miRNA-pathway relationship network construction

Construction of ceRNA network

Discussion

Conclusion

Methods

Tissue Samples and RNA Extraction

Next generation sequencing and data processing

The association analysis between lncRNA and mRNA

MiRNA target prediction and Construction of ceRNA network

Functional enrichment and Visualization of ceRNA network

Abbreviations

Declarations

Acknowledgements

Funding

Availability of data and materials

Authors’ contributions

Ethics approval and consent to participate

Consent for Publication

Competing interests

Reference

Tables

Supplementary Files

Status:

Version 1