LncRNAs Function as Critical Metabolic Regulators in yak Liver by Integrative Transcriptome Analysis

Background: The yak (Bos grunniens) is a crucial resource to supply meat and milk to the people localized in Qinghai-Tibetan plateau area. To identify lncRNAs regulating metabolism in yak, this work adopted transcriptome method to simultaneously prole mRNAs and lncRNAs of liver in yak under three representative age (LD: Liver 1 Day, LM: Liver 15 Months, LY: Liver 5 Years) conditions. Result: Of 288 differentially expressed lncRNAs, function-oriented selection yield 88 regulated metabolically related lncRNAs that were differentially expressed at least two age conditions. These lncRNAs predicted by lncRNA-mRNA correlation analysis to function in various aspects of metabolism. Selected regulations of liver metabolically related lncRNAs were further veried by qRT-PCR. Conclusion: Combining high throughput RNA-seq screening screens, bioinformatics predictions, lncRNA-mRNA correlation analysis and qRT-PCR analysis, this study supports that a class of lncRNAs function as important metabolic regulators and establishes a foundation for further investigating the role of lncRNAs in yak.


Background
The yak (Bos grunniens) lives throughout the Qinghai-Tibetan plateau in western China at high altitudes where harsh and variable plateau climate exist including low humidity, temperature, and oxygen levels; strong winds and ultraviolet radiation. The yak has adapted to thrive under these harsh conditions, where few other livestock could survive. The yak is a crucial resource to supply meat and milk to the people localized in Qinghai-Tibetan plateau area. For these reasons, the yak is one of the most critical domestic animals for the 6.5 million Tibetans [1]. Based on the fact that the yak represents the primary source of meat and milk for the people of this region, it is imperative to understand the genetic makeup of the yak and how this can in uence its development.
While, metabolic process is one of the most basic activity which is related to all fundamental biological and physiological activities in all animals [2]. Liver, the most important metabolic organs, engage in continuous dialogs to regulate overall metabolic balance by endocrine factors. And it is also a complex digestive gland in ruminants including yak and plays a critical role in the metabolism of substances.
More importantly, during the developmental period, liver plays a critical role [3][4][5]. Dorland et al. previously reported the period of transition from new born to young cattle involved considerable metabolic adaptation in dairy cows and the liver [3]. Therefore, in the present study, the liver of new born (1 day), young yak (15 months) and mature yak (5 years) were chosen to be sampled.
Although a lot of work has been done in studying the regulation of individual metabolic pathways in recent years, there are still quite a few unknow elds in the complicated regulatory networks that control metabolic physiology. It remains di cult to develop useful treatment strategy against metabolic disorders such as diabetes and obesity [6][7][8]. Nearly two-thirds of transcripts are noncoding RNAs, mainly from regions previously considered as junk gene [9]. Among the noncoding RNAs identi ed so far, lncRNA is the longest and most di cult to understand, which are transcripts of 200 nt or longer that lack protein coding potential. LncRNAs have been recognized in all model organisms [10], and their number has been increasing continuously [11][12][13]. It has been shown that lncRNAs have an in uence on several aspects of cellular function, including chromatin modi cation and transcription regulation, RNA stability, and translational control [14,15]. Evolution studies have shown that more than 1,000 lncRNAs may have conservative functions in mammals [16]. There were also many reports in bovine indicates lncRNAs function in embryo development [17], skeletal muscle [18] and mammary glands [19]. In mouse, in silico method had been used to identify functional lncRNAs in metabolic process [20]. However, it is still unknown elds in ruminant animals like bovine or yak. Therefore, in order to systemically illustrate the importance of lncRNAs in the homeostasis of metabolism of yak, this work combined a high throughput RNA-seq screening of lncRNAs transcripts from liver of yak in diverse ages to establish a comprehensive work ow to determine functional lncRNAs in metabolic regulation. LncRNA-mRNA network analysis is used to predict lncRNAs function, and qRT-PCR analysis connects lncRNAs with speci c metabolic pathways. In conclusion, our research supports lncRNAs as a critical composition of metabolism and provides a fundation for systemically identifying and predicting functional lncRNAs that regulate the homeostasis of metabolism.

Result
Dynamic regulation of lncRNAs and mRNAs in liver under different age conditions Both mRNA and lncRNA transcripts were detected in our study, current expression pro le could compare the transcripts of lncRNAs and mRNAs at different developmental age stages in the yak liver. Hierarchical cluster analysis was performed on all expressed transcripts, which showed the mRNAs expression pro les clearly divide all samples into three different groups according to age stages (1 day, 15 months, 5 years) and these samples were closely clustered together in each group of different developmental age stages (Fig. 1A, top left). But what is interesting to us is that an almost identical sample aggregation pattern showed in the result of lncRNAs (Fig. 1A, top right), which indicates that expression pro les of lncRNAs could be used as a signal indicating similar to those of protein-encoded mRNAs. Using the RNAseq method, 35216 mRNAs and 10073 lncRNAs transcripts were detected from transcripts in the yak livers at three different developmental age stages (LD: Liver 1 Day, LM: Liver 15 Months, LY: Liver 5 Years) ( Fig. 2B).
In addition, all samples were divided into diverse groups for both mRNAs and lncRNAs through PCAs analysis on all regulated transcripts (Fig. 1C), which indicate that regulated lncRNAs and mRNAs transcripts might coordinate related biological processes. Furthermore, to determine their functional connections, analysis of the network with lncRNA-mRNA through related analysis of samples were conducted. 433 mRNAs and 152 lncRNAs were identi ed to be regulated by LD versus LM, 412 mRNAs and 160 lncRNAs by LD versus LY, as well as 263 mRNAs and 66 lncRNAs by LM versus LY (Fig. 1D), suggesting that their expression levels were strictly controlled by age conditions of the yak. Functional analysis of differentially expressed genes (DEGs) in each age condition of yak liver At LD versus LM, the DEGs of many functional Gene ontology (GO) categories enriched for metabolism, ion binding, developmental process and so on. While in KEGG pathway analysis, DEGs were also enriched for metabolism, biosynthesis, ECM-receptor interaction and so on ( Fig. 2A). At LD versus LY, there were functional GO categories enriched for metabolism, tissue remodeling and developmental process and so on. While in the KEGG pathway analysis, DEGs were also enriched for terms like metabolism, biosynthesis, PI3K-Akt signaling pathway and so on (Fig. 2B). Similar functional GO categories of metabolism, biosynthetic process and oxidation reduction process and so on were enriched at LM versus LY. While in the KEGG pathway analysis, DEGs were also enriched for terms like metabolism, biosynthesis, focal adhesion and so on (Fig. 2C). All these results suggest metabolism related function were developed during the aged process of yak (Fig. 2).
Functional prediction of metabolically related differentially expressed lncRNAs by lncRNA-mRNA co-expression correlation under at least two age conditions By correlation of lncRNA-mRNA expression interaction network, dynamically regulated lncRNAs were predicted their function at different developmental ages. 88 non-redundant lncRNAs regulated at least two of the three age conditions were selected from the totally 288 regulated differentially expressed lncRNAs (Fig. 1D). Furthermore, to detect the speci c function of these metabolically associated lncRNAs as example, six metabolically related lncRNAs were selected for further analysis (Fig. 3).
Taking six lncRNAs as examples to make the functional analysis, correlation analysis predicts their role in metabolism or cell differentiation. ENSBMUG00000000490 and XLOC_045379 were predicted to be related with fat single organism metabolic process ( Fig. 3A and 3B), and XLOC_021536 and XLOC_041441 associated with collagen metabolic processes and protein metabolic processes in liver, respectively ( Fig. 3C and 3D). Moreover, ENSBMUG00000026019 and XLOC_183608 have been shown to be enriched in cell proliferation and transport ( Fig. 3E and 3F). These results indicate that the established co-expression network can effectively predict the potential metabolic functions of age-regulated and metabolically sensitive lncRNAs. In addition, the dynamic regulation of several randomly selected lncRNAs under different age conditions was con rmed by qRT-PCR, proving the feasibility of this method (Fig. 4).

Discussion
Integrative lncRNAs function related to metabolic were selected using high-throughput RNA-seq screening, and a work ow for the discovery and characterization of functional lncRNAs were established in yak at different ages. The data obtained in this experiment support that lncRNAs was a crucial component of metabolism. Both the expression pro le of lncRNAs and RNAs have change coordinately according to different developmental ages including 1 day,15 months, and 5 years. Moreover, groups of metabolism relative lncRNAs were regulated in liver by ages, and their expression often changes signi cantly in yaks at diverse ages, supporting their potential biological signi cance.
Identifying lncRNAs regulating metabolic and nding its functional properties is still a hard task in animals and only a few lncRNAs could regulate metabolism were reported for now [21][22][23][24]. Moreover, it is impossible to directly determine the role of lncRNAs in metabolic process, for the reason of its non-coding feature. In order to get over this di culty, this article build a functional lncRNAs detection method by interactive analysis of lncRNAs in liver of yaks at different ages. This integrated approach e ciently reduce 288 age-adjusted lncRNAs to 88 putative metabolic lncRNAs which is differentially expressed at least in two age stages.
Only a few previous published papers showed that lncRNAs play a critical role in regulating metabolic pathway in mice [21][22][23][24], and certainly there is no report on the function of lncRNAs in yak. In order to understand the processes of yak liver metabolism in depth, more lncRNAs need to be discovered and characterized, and con rming the work ow of lncRNAs could speed up the selection and characterization of metabolic lncRNAs, which offered more perspective for the complex metabolism networks.
Making out the impact of lncRNAs in metabolic pathways at different developmental ages could enhance the unknown elds of complex metabolic physiology. Comparative genomics method had been widely used in revealing potential functions of novel protein-coding genes based on information of homologous genes or protein motifs, but has proved to be noneffective in nding the function of lncRNAs [25]. For most lncRNAs have not been functionally studied, and more importantly, lncRNAs were much more conservative than mRNAs, even between closely related model organisms [26]. Since most lncRNAs are unique to different animals, it is di cult to infer their role according to the sequence matching or evolutionary records [27].
For it is hard to predict functions for lncRNAs, it becomes very di cult to discover and characterize lncRNAs that regulate metabolism. Since metabolism is essential for almost all the biological processes, all organisms, any controlling with vital metabolic pathways is achieved usually through reconnecting metabolic uxes. Therefore, a strict functional measurement is needed to determine critical points in metabolic regulation. To check the function of lncRNA further in animal models, it is rstly necessary to correctly infer the functional information of lncRNAs. In addition, the results of this experiment can also be bene cial to design targeted and detection methods to determine the unique metabolism in speci c lncRNA, which is usually suppressed in complex interaction or compensation situation. Therefore, an e cient method to predict metabolic lncRNAs or connecting lncRNAs to metabolic pathways could signi cantly accelerate the identi cation of important lncRNAs metabolic regulators.

Conclusions
The present study has illustrated the expression of metabolism-related lncRNAs in different ages of yak and lay a foundation for future nding and characterizing metabolism relative lncRNAs. Potential speci c lncRNAs metabolic regulators with possible functions has also been predicted. This study could provide useful information for revealing the metabolic functions of lncRNAs in yak liver.

Animal
Totally nine Maiwa yaks from Yak Technology Park (Hongyuan County of Sichuan Province in China) at diverse age conditions (1 day, 15 months and 5 years old) were used in this study and each developmental age stage had three yaks. The experimental animals were healthy and under the same management. Yaks were weighing 10.55 ~ 14.24 kg, 96.38 ~ 101.37 kg and 240.73 ~ 296.36 kg. All yaks were stunned with a captive bolt pistol (Cash 8000 Model Stunner, 0.22 calibre, 4.5 grain cartridge) to ameliorate the suffering of the animals prior to their humane killing, following which exsanguination via a transverse incision of the neck was carried out in the slaughterhouse. Then, the liver tissues were excised immediately and rapidly stored in liquid nitrogen until RNA isolation. Yaks with different age were established for pro ling the transcriptomes in liver tissue.

RNA quanti cation and quali cation
RNAiso Plus (TaKaRa, Japan) was used to extract total liver tissue RNA. The extent of RNA degradation and contamination was monitored on 1% agarose gels. RNA purity was checked using the NanoPhotometer ® spectrophotometer (IMPLEN, CA, USA). RNA concentration was measured using Qubit ® RNA Assay Kit in Qubit ® 2.0 Flurometer (Life Technologies, CA, USA). RNA integrity was assessed using the RNA Nano 6000 Assay Kit of the Bioanalyzer 2100 system (Agilent Technologies, CA, USA). If the OD 260 / 280 value of total liver tissue RNA is between 1.8 ~ 2.2, the quality inspection is quali ed Table 1) and the next analysis can be carried out.
Library preparation for high-throughput sequencing A total amount of 3 μg RNA per sample was used as input material for the RNA sample preparations. A chain-speci c library was constructed by removing rRNA, that was, sequencing libraries were generated using NEBNext® UltraTM RNA Library Prep Kit for Illumina® (NEB, USA) [28]. Brie y, mRNAs were puri ed from total RNA using poly-T oligo-attached magnetic beads. First strand cDNA was synthesized using random hexamer primer and RNase H. Second strand cDNA synthesis was subsequently performed using DNA Polymerase I and RNase H. The puri ed double-stranded cDNA was subjected to end repair, was added poly-A tail, and was connected to a sequencing adapter. In order to select cDNA fragments of preferentially 150~200 bp in length, the library fragments were puri ed with AMPure XP system (Beckman Coulter, Beverly, USA). Then 3μL USER enzymes (NEB, USA) were used to degrade the second strand of U-containing cDNA, and perform PCR ampli cation to obtain the library. After the library was constructed, Qubit2.0 Fluorometer (Invitrogen, USA) was used for preliminary quanti cation. PCR products were puri ed (AMPure XP system) and library quality was assessed on the Agilent Bioanalyzer 2100 system. qRT-PCR (Bio-Rad, USA) was used to accurately quantify the effective concentration of the library to ensure library quality. The library preparations were sequenced on an Illumina Hiseq platform and 150 bp paired-end reads were generated.

Transcriptome sequencing data analysis
The quality-controlled Clean Reads was compared with the yak reference genome (BosGru_v2.0: ftp://ftp.ensembl.org/pub/release-97/fasta/bos_mutus/) quickly and accurately by the HISAT2 v2.0.5 software (http://ccb.jhu.edu/software/hisat2) to obtain Mapping Reads for subsequent analysis [29]. At the same time, the comparison results were evaluated for quality. By analyzing the different regions and chromosome distributions of Mapping Reads in the reference genome, to obtain the comparison e ciency and mapping information about Mapping Reads for each sample [30]. The new transcripts for Mapping Reads were assembled and quanti ed by StringTie (1.3.3b) software [31]. Sequencing depth and gene length were corrected using FPKM [32]. Differential expression analysis was performed using the DESeq2 R package (1.10.1). DESeq2 provide statistical routines for determining differential expression in digital gene expression data using a model based on the negative binomial distribution. Genes with an adjusted padj<0.05 found by DESeq2 were assigned as differentially expressed [33].

Functional enrichment
ClusterPro ler (3.4.4) software was used to perform GO function enrichment analysis and KEGG pathway enrichment analysis on differential gene sets, to predict the biological processes and functions that they may participate in, and to make corresponding classi cation and statistics [34]. All lncRNAs were used for target gene prediction, that was, the target genes of lncRNAs were predicted by the positional relationship (co-location) and expression correlation (co-expression) of lncRNAs and protein-coding genes [35]. Then, functional enrichment analysis (GO / KEGG) was performed on the target genes of differential lncRNAs to predict the main function of lncRNAs.
Real-time uorescence quantitative PCR Measurement of gene expression with qRT-PCR has been applied in our studies. Brie y, total RNA was extracted from liver tissue of each group using a RNAiso Plus (TaKaRa, Japan). RNA was complete without degradation, good quality and high purity, which meet the requirements of subsequent experiments. qRT-PCR was performed in a total volume of 10 µL containing 5.2 µL of TB Green TM Premix Ex Taq TM (TaKaRa, Japan), 1 µL of cDNA, 0.8 μL of each primer (10 mM) and 2.2 µL of double-distilled water. The reaction condition used were as follows: 95 °C for 3 min, followed by 39 cycles of 95 °C for 10 s, 53.4/52 °C (β-actin/veri cation genes) for 20 s and 72 °C for 20 s, with the dissolution curve increasing from 0.5 °C to 95 °C every 5 s. The assays were performed on a real-time uorescence quantitative PCR System (Bio-Rad, USA). For each sample, the cycle threshold (CT) values were obtained from three replicates. β-actin mRNA was employed as an internal reference. The primers used for ampli cation of target and internal reference genes were presented in Supplementary Table 2

Consent for publication
All of the authors have approved the nal version of the manuscript, agree with this submission to BMC Genomics.

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