Detection and integrated analysis of lncRNA and mRNA relevant to plateau adaptation of Yak.

Known as the 'ship of the plateau', through thousands of years evolution and cruelty environments selection containing low oxygen and strong ultraviolet radiation, yaks have adapted plateau environments and supplied important goods and materials for the people in the Qinghai-Tibet Plateau. This study aimed to identify differentially expressed (DE) genes and novel long non-coding RNAs (lncRNAs) of yaks for the Plateau adaptation and their underlying co-expression and regulatory network. We carried out RNA-seq analysis for cerebral and cerebellar tissue specimens of Bos taurus, Bos grunniens × Bos Taurus and B. grunniens. Furthermore, 12,072 pseudo lncRNAs were predicted using three software. In total, 4,257 significant DE transcripts were identified using the Ballgown R package (p < .01), of which 1,021 were protein-coding genes, 14 were known lncRNAs, and 661 were novel lncRNAs. Using WGCNA, a co-expression network of DE mRNAs and lncRNAs comprising 5 modules was generated to determine functional associations clusters. This study reveals a valuable sub-network comprising 8 hub genes, one known lncRNA and 5 novel lncRNAs in the major module. These hub genes are associated with blood pressure regulation, generation of reactive oxygen species and metabolism. The analysis of co-expressed genes thus provides a basis for the regulatory mechanisms in PA in Yaks and for the detection of additional genes between cross-breed and parent populations.

INTRODUCTION

Long non‐coding RNAs (lncRNAs; >200 nt in length) regulate translation or protein synthesis (Beck et al., 2014; Cloutier, Wang, Ma, Petell, & Tran, 2013; Cusanelli & Chart rand, 2015) as well as mRNA translation in different environments (Riaz, Wolden, Gelblum, & Eric, 2016). The lncRNA‐mRNA network regulates environmental adaptation and biological responses (Huang et al., 2014; Salleh, Mazzoni, Løvendahl, & Kadarmideen, 2018). To further elucidate this network, the general linear model has been used to analyse the underlying associations.

In the Qinghai‐Tibet Plateau, the Yak (Bos grunniens) is a prominent livestock species to the Tibet region (Hu et al., 2012; Qiu et al., 2015). Yaks are called the ‘ship of the plateau’ and are highly adaptable to low the atmospheric pressure, hypoxia and high ultraviolet radiation in the plateau (Hu et al., 2012; Jincheng et al., 2014; Zhang et al., 2016). To assess the genetic advantages of yaks, numerous studies have reported associations of SNPs and copy numbers variations with plateau adaptation (PA; Hu et al., 2012; Medugorac et al., 2017; Qiu et al., 2015; Zhang et al., 2019); however, the underlying gene expression, regulation and the functions of PA remain unclear, especially the mRNA‐lncRNA co‐expression network. The brain is a key tissue as a nerve centre and is responsible for physiological response to adaptation in extreme environments (Marshall‐Goebel, Damani, & Bershad, 2019). Most of the cattle (Bos Taurus) died from intracranial haemorrhage after exposure to high altitudes of over 3,000 m, but the Yak survived very well and were highly useful to the people in the plateau (Gerald, Han, & Long, 2003). These evidences indicated the brain is an important and complex system for mammal adaptation of plateau environments. In the 2016, Wang et al. used four organs of Yak and Cattle to provide insights into high‐altitude adaptation (Wang et al., 2016). But because of the large gap and differences in genetics between B. grunniens and B. taurus, whatever chosen any one reference genomes the noise of population was always unable to distinguish. Actually, as the cross‐breed of Yak (B. grunniens) and Cattle (B. taurus), Zang Yakow presents moderate PA. With the help of expression trend in three breeds, it will improve the accuracy to detect candidate genes and predict the co‐expression network involved in PA (Ma et al., 2015; Wu et al., 2018).

In this study, we performed RNA‐seq analysis for three breeds: Leiwuqi Yak (B. grunniens), Zang Yakow (B. grunniens × B. taurus) and Sanjiang Cattle (B. taurus). The differences in the PA potentials among these three breeds were determined through high‐throughput transcriptome sequencing of cerebral and cerebellar tissues from animals to identify DE genes and predict novel lncRNAs. Furthermore, we constructed a network of DE mRNAs and lncRNAs to assess their interactions and the role of these interactions.

MATERIALS AND METHODS

Sample collection and RNA extraction

All animals had a similar health status and were approximately 4.5 years old and females. Eighteen samples (the cerebrum and cerebellum came from 3 Leiwuqi Yak, 3 Zang Yakow and 3 Sanjiang Cattle) were harvested in the Leiwuqi region (Leiwuqi Yak and Zang Yakow) and Wenchuan region (Sanjiang Cattle), Sichuan, China. The brain cortex of all samples was collected after slaughter while considering animal welfare in October 2018. Three individuals of the same breed were randomly selected as replicates in a herd of >50 individuals. The Leiwuqi Yak was considered the M group (living average 4,450 m altitude); Sanjiang Cattle was considered H group (living average 1,325 m altitude); Zang Yakow was considered Z group (living average 3,200 m altitude). All clean tissue samples were divided into 1–2 cm3 sections and snap‐frozen in liquid nitrogen. Using TRIzol® Plus RNA Purification Kit (Invitrogen) in accordance with the manufacturer's instructions, total RNA was extracted from brain tissue. Ribosomal RNA (rRNA) was eliminated before sequencing, using the Epicentre Ribo‐zeroTM rRNA Removal Kit in accordance with the manufacturer's instructions.

RNA sequencing and transcript assembly

After eliminating rRNA, cDNA sequencing libraries of total coding RNAs and ncRNAs were generated using the mRNA‐seq sample preparation Kit (Illumina). Using the Illumina HiSeqTM 4000 platform, we obtained paired‐end sequencing reads of 18 tissue samples. The sequencing reads were filtered using fastp (version 0.19.8) software with default parameters. Clean reads were aligned to the Yak reference genome (GCF_000298355.1 BosGru v2.0; Qiu et al., 2012) with HISAT2 (version 2.1.0) to generate a Sam file (Kim, Langmead, & Salzberg, 2015). After sorting and conversion to bam with Samtools (version 1.9), we used StringTie (version 1.2.3; Pertea et al., 2015) to assemble transcripts for each sample. All transcripts were merged with the ‘StringTie‐merge’ function of StringTie with the corresponding sample file names, and output fields were generated in Ballgown format (Frazee et al., 2015). The merged file was annotated with the reference genome annotation file (GCF_000298355.1 BosGru v2.0) via gffcompare (version 0.10.8).

Identification of lncRNAs

After pre‐processing the RNA‐seq data, all transcripts in the merged transcript annotation file were identified using four filters to predict potential transcript. First, transcripts ≤200 nt were eliminated. Second, based on the relative position among transcript and known genes, we used gffcompare to mark the transcript type. Those transcripts with class code ‘u’ (unknown intergenic), ‘x’ (exonic overlap on the opposite strand), ‘i’ (fully contained within a reference intron), ‘j’ (multi‐exon with at least one junction match) and ‘o’ (other same strand overlap with reference exons) were retained. Third, transcripts with low expression levels (FPKM ≤ 1) were filtered out. Fourth, transcripts more than 1 exon were retained (Cloutier et al., 2013; Fan et al., 2019). The FASTA format data were extracted using the getfasta function of bedtools from the reference genome (Quinlan & Hall, 2010). Whole remaining transcripts were considered new transcripts. Thereafter, three methods were used to predict the potential protein‐coding potential of these new transcripts. These three methods are the Coding Potential Calculator (CPC2; Kang et al., 2017), the Coding‐Non‐Coding‐Index (CNCI; Sun et al., 2013) and the Coding Potential Assessing Tools (CPAT; Wang et al., 2013). The intersection of these three predicted lncRNAs was considered to comprise the new lncRNAs for downstream analyses.

DE mRNAs and lncRNAs among M and H groups

New lncRNAs, known mRNAs and known lncRNAs were identified from the reference genome in accordance with the annotation file and through identification of new lncRNAs. The Fragment PerKilobase Million (FPKM) metric was used as the unit of measurement to estimate transcript abundance, using R package Ballgown, which is an integrated tool for statistical analysis of assembled transcriptomes. The group was considered a covariate factor, and different tissues (cerebellum and cerebrum) were considered to adjust the variate factors in the Ballgown model for DE analysis. The thresholds |log2FoldChange|>1 and p < .01 were assigned for DE genes, using the Ballgown package.

GO enrichment analysis was conducted to cluster genes in accordance with the biological process, cellular component and molecular function. Candidate DE genes were analysed using R package “clusterProfiler” with database “org.Bt.eg.db” (R Development Core Team, 2016; Yu, Wang, Han, & He, 2012) for statistical analysis and visualization of functional profiles of genes and gene clusters.

Co‐expression network for DE RNA and PA

The expression levels of these significant DE genes among the whole three populations (Yak, Yakow, and Sanjiang Cattle) were determined to construct a co‐expression network using the weighted gene co‐expression network analysis (WGCNA) R package (Langfelder & Horvath, 2008). The WGCNA method is a systems biology method for describing the correlation patterns among genes across microarray samples. It can be used for clusters (modules) of highly correlated genes, for summarizing such clusters using the module eigengene or an intramodular hub gene, for relating modules to one another and to external sample traits (using eigengene network methodology), and for calculating module membership measures (Langfelder & Horvath, 2008). The soft power was set from 1 to 30 lops to optimize the association model. Among all genes and lncRNAs, pairwise Pearson's correlations were determined to generate an adjacency matrix. The topological overlap measure (TOM) was estimated via the optimum model using the adjacency matrix.

A heat map was generated with these module eigengenes, using the first principal component to determine variations among genes and lncRNAs expression. All genes in each module were considered to visualize network associations with Cytoscape (Shannon et al., 2003). Major hub genes in the largest module were selected to enrich with correlated known lncRNAs and novel lncRNAs.

qPCR validation

RNA samples, which were extracted from the same sources and at the same concentration as library preparation, were reverse‐transcribed to cDNA using the PrimeScript™ RT reagent Kit with gDNA Eraser (Perfect Real Time; TaKaRa, Dalian China) according to the manufacturer's instructions. Four hub genes in the most related module of the co‐expression network were selected to validate RNA‐seq result. Primers were designed using Primeprimer 5.0 (Table S2). As a reference gene, cattle GAPDH was used to adjust relative CT value. Relative quantifications of genes were calculated by the ∆CT method (Rao, Huang, Zhou, & Lin, 2013).

RESULTS

Sample collection and RNA sequencing

Each sample was totally mapped at an average mapping ratio of 91.57% relative to the reference genome. Of these, 84.22% were uniquely mapped reads and the remaining 2.98% were multiple mapped reads (Table S1).

A summary pipeline of the experimental workflow, bioinformatics and statistical analysis is displayed (Figure S1). In total, 101,835 transcripts were assembled in the merged assembled file of 18 samples, of which there were 32,667 transcripts expressed in all 18 samples. All transcript expression levels in the 18 samples were presented (Figure S2) with log2 (FPKM + 1) values. Furthermore, we determined the distribution of transcripts per known gene (Figure S3). There were 3,755 genes with only one transcript, and the ratio of total genes was 24.16%, indicating that most genes had multiple transcripts. We designed a smoothened colour density representation (Figure S4) to verify the repeatability and compare the similarities and differences in the expression levels between the cerebrum and cerebellum in each sample.

55,542 new transcripts were obtained to predict the protein‐coding potential with three software. CPC2 predicted 21,642 pseudo lncRNAs; CNCI predicted 16,819 pseudo lncRNAs; CPAT predicted 36,759 pseudo lncRNAs. Of these, 12,072 pseudo lncRNAs were detected at the intersection of the aforementioned three software (Figure 1a) and were considered novel lncRNAs. Most of the novel lncRNAs were unknown intergenic (64.17%), 19.75% were multi‐exon matched, 10.32% were fully contained within introns, 3.98% overlapped with opposite sites, and 1.78% overlapped with exons (Figure 1b). Comparisons between the FPKM distribution of novel lncRNAs and known lncRNAs are presented in Figure 1c.

FIGURE 1

Genomic features and classification of novel lncRNA. (a) The intersection of predictive long non‐coding RNA by three methods (CNCI, CPC2 and CPAT). (b) Classification of novel lncRNA according to ‘class code’ showing the type of relative position between a transcript and the closest reference transcript. (c) The FPKM distribution of Novel lncRNA and known lncRNA

Analysis of DE genes

To identify DE genes between Yak (M) and Sanjiang Cattle (H), we filtered out those transcripts of no variational expression in the linear model. In total, 4,257 significant (p < .01) genes were identified using Ballgown R package with species as covariate factors and tissues as adjustable variate factors. After removed |log2foldchange|<1, there were 2,577 different expression (DE) genes, and 867 of DE genes were upregulated and 1,710 were downregulated (Figure 2a). Of these, 1,021 were protein‐coding genes, 14 were known lncRNAs, and 661 were novel lncRNAs. The DE genes were identified in accordance with the top 5 GO terms (Figure 2b). There was few DE genes between the cerebrum and cerebellum in all three groups (Figure S5) that indicated the differential express of these two tissues. We also compared the DE genes between the Yakow and Sanjiang Cattle (Figure S6). There were 8,517 DE genes.

FIGURE 2

The screening and enrichment of Different Expression genes in the Yak and Cattle group. The DE genes were filtered between Yak and Cattle with 2 times fold change (a). The 14 GO terms of all subontologies were significant enriched with whole DE genes (b)

Co‐expression network

The R square of the filter criteria was set to 0.8 to optimize and select the soft power of the prediction model (Figure S7). We also used the scaleFreePLot function of WGCNA to verify the relative levels between the whole estimated network and scale‐free network (Figure S8 and Supplementary file lnc_known_mRNA_node.txt). The associated expression genes were clustered into 6 groups (the grey module was the mismatched module), and each group (excluding the grey module) contained more than 50 genes (Figure 3). Based on the TOM, the correlation of genes in the all modules was visualized (Figure 4; Figure S9) via Cytoscape. The most of genes in the largest module (turquoise) were enriched into 6 GO terms (Figure S10). Comparison with previous study, there are four important pathways related to PA (Qi et al., 2019; Song et al., 2019; Xiong et al., 2015). We used KEGG pathway database to match our DEGs to these four pathways (Figure S11).

FIGURE 3

Optimum modules with the whole coding genes, known lncRNA and novel lncRNA. (a) The cluster dendrogram of whole DE genes. The total of five modules was created with the expression of coding gene, known lncRNA and novel lncRNA. (b) The relationship and among five modules

FIGURE 4

Co‐expression network in modules with the whole coding gene, known lncRNA and novel lncRNA. (a) Co‐expression network visualization: nodes show coding gene (red), known lncRNA (blue) and novel lncRNA (green). (b) Classification of RNA in the whole co‐expression network

In total, 979 edges (genes and lncRNAs) were associated with the most significant terms of the largest related module. We selected the top 10 hub genes, 5 known lncRNAs and 5 novel lncRNAs to elucidate the underlying regulatory functional network based on the hub genes (Figure 5). All lncRNAs and mRNAs expression of whole network were used to calculate Neighbour‐Joining tree (NJtree) and principle components analysis (PCA) to present population structure (Figures S12 and S13). And the variance explained proportion of top two PCs was more than 70% (Figure S14).

FIGURE 5

Sub‐network with most enrichment genes and lncRNAs. The top 10 enrichment hub genes (Blue Diamond) in whole network, top 5 known lncRNAs (Green Circle) and top 5 novel lncRNAs (Red Triangle) were selected to build sub‐network. The more enrichment line indicated more related regulation of expression

The related expression levels of random selected 4 hub genes (MRPL49, ERAP1, IL34 and PINK1) were verified by quantitative reverse transcription PCR (RT‐qPCR). The correlations between relative CT value and RNA‐seq expression of 4 hub genes were 0.67 ~ 0.99 in cerebrum (Figure 6) and 0.47 ~ 0.97 in cerebellum (Figure S15).

FIGURE 6

RT‐qPCR validation with four hub genes expression in Cerebrum. The RT‐qPCR result of MRPL49 (a), IL34 (b), ERAP1 (c) and PINK1 (d) was used to verify RNA‐seq expression result

DISCUSSION

Genetic characteristics of Yak, Yakow and Cattle

Yak and Cattle were estimated to have diverged approximately 4.9 million years ago; therefore, the Yakow cross‐breed between Yak and Cattle presented a male sterile phenotype, which was gradually restored with a back‐cross between female Yakow and Yak after several generations (Xi et al., 2012). In the 2016, Wang et al. used four organs of Yak and Cattle to provide insights into high‐altitude adaptation (Wang et al., 2016). This manuscript is the first time to detect Plateau Adaptation genes from Yak, Yakow and Cattle population. The reads unique mapping ratio (Table S1) of Yak was the highest, followed by that were Yakow and Cattle, thus indicating the associations of their genetic backgrounds relationship with the reference genome.

In the WGCNA, the clustered genes in one module had the same expression trend among all groups. The 10 hub genes in turquoise module were enriched in macrophage colony‐stimulating factor receptor binding and atom of oxygen terms. Of these hub genes, ERAP1 is a key gene significantly contributing to the regulation of blood pressure in mammals (Hisatsune et al., 2015). Oxygen stress in the endoplasmic reticulum is decreased upon a reduction in ERAP1 levels, and blood pressure is decreased to maintain a balance (Goto, Ogawa, Hattori, & Tsujimoto, 2011). IL34, GRAMP3, PINK1 and MRPL49 are reportedly associated with the regulation of the generation of reactive oxygen species and metabolism (Geisler et al., 2010; Jung, Seo, Lee, Lee, & Roe, 2015; Radaeva & Simbirtsev, 2019; Wang et al., 2012; Xiao et al., 2017).

It was reported in 2012 that keratinocytes and neurons are the main sources of IL34 and that IL34 specifically directs the differentiation of myeloid cells in the skin epidermis and central nervous system (Wang et al., 2012). It is both a cytokines and a vasoactive peptides, and high expression levels can reduce the rate of complications (myocardial infarction, acute cerebrovascular events and transient ischaemic attacks; Radaeva & Simbirtsev, 2019). IL34 was highly expressed in the Yak group but its expression was very low in the Cattle group according to RNA‐seq and qPCR results. In 2018, Lan sequenced Jinchuan Yak samples and found that PINK1 was enriched in the Parkinson's disease pathway (Lan et al., 2018). The PINK1 expression levels of Yak were significantly higher than those of Yakow and Cattle. Sarkar Soma used gene expression to reveal that MRPL49 is associated with the hypoxia due to high‐altitude acclimatization (Soma, 2012). In 2017, Mota and Guimaraes reported that MICAL1 affected age at the first calving in the Nellore cattle (Braeutigam et al., 2014). Reestablishment of elevated TLE1 levels by ectopic expression protects neurons from death, whereas the suppression of TLE1 expression in otherwise healthy neurons induces cell death (Dastidar, Narayanan, Stifani, & D’Mello, 2012). This maybe the key reason for the intracranial haemorrhage observed after exposing Sanjiang cattle to a high altitude of over 3,500 m.

CONCLUSION

The present results indicate that the PA is regulated by numerous genes and lncRNAs, which comprise an underlying network. The functions of these hub genes are associated with the regulation of blood pressure, generation of reactive oxygen species and metabolism. The analysis of co‐expressed genes with the PA provides a basis to understand the regulatory mechanisms in the complex PA of Yak.

COMPETING INTERESTS

The authors declare that they have no competing interests.

AUTHORSHIP STATEMENT

JBW: Processed and analyzed RNA sequencing data, implemented software, method and drafted the manuscript; LD: Processed qPCR; ZXC: Participated in discussions regarding data analyses; JKW: Participated in discussions regarding interpretation of results; HW: Performed animal phenotyping and analyses for animal selection, collected samples; YT: Participated in discussions regarding interpretation of results; JCZ: Designed experiments, participated in discussions regarding data analyses. All authors made significant contributions editing the manuscript. All authors read and approved the final manuscript.

ETHICS APPROVAL AND CONSENT TO PARTICIPATE

The protocol for Animal Care and Use was approved by the Animal Ethics and Welfare Association of the Southwest Minzu University (No. 16053), and experiments were performed according to the regulations and guidelines established by this committee.

Supplementary Material

Click here for additional data file.