Identification, functional prediction, and key lncRNA verification of cold stress-related lncRNAs in rats liver

Abstract

Cold stimulation reduces the quality of animal products and increases animal mortality, causing huge losses to the livestock industry in cold regions. Long non-coding RNAs (lncRNAs) take part in many biological processes through transcriptional regulation, intracellular material transport, and chromosome remodeling. Although cold stress-related lncRNAs have been reported in plants, no research is available on the characteristic and functional analysis of lncRNAs after cold stress in rats. Here, we built a cold stress animal model firstly. Six SPF male Wistar rats were randomly divided to the acute cold stress group (4 °C, 12 h) and the normal group (24 °C, 12 h). lncRNA libraries were constructed by high-throughput sequencing (HTS) using rat livers. 2,120 new lncRNAs and 273 differentially expressed (DE) lncRNAs were identified in low temperature environments. The target genes of DElncRNA were predicted by cis and trans, and then functional and pathway analysis were performed to them. GO and KEGG analysis revealed that lncRNA targets were mainly participated in the regulation of nucleic acid binding, cold stimulation reaction, metabolic process, immune system processes, PI3K-Akt signaling pathway and pathways in cancer. Next, a interaction network between lncRNA and its targets was constructed. To further reveal the mechanism of cold stress, DElncRNA and DEmRNA were extracted to reconstruct a co-expression sub-network. We found the key lncRNA MSTRG.80946.2 in sub-network. Functional analysis of key lncRNA targets showed that targets were significantly enriched in fatty acid metabolism, the PI3K-Akt signaling pathway and pathways in cancer under cold stress. qRT-PCR confirmed the sequencing results. Finally, hub lncRNA MSTRG.80946.2 was characterized, and verified its relationship with related mRNAs by antisense oligonucleotide (ASO) interference and qRT-PCR. Results confirmed the accuracy of our analysis. To sum up, our work was the first to perform detailed characterization and functional analysis of cold stress-related lncRNAs in rats liver. lncRNAs played crucial roles in energy metabolism, growth and development, immunity and reproductive performance in cold stressed rats. The MSTRG.80946.2 was verified by network and experiments to be a key functional lncRNA under cold stress, regulating ACP1, TSPY1 and Tsn.

Figure 1

Informatics pipeline for the identification of lncRNAs.

Figure 2

Genomic features of lncRNAs. (A) Raw data composition. (B) Distribution of lncRNAs on each chromosome. (C) Venn diagram to predict the coding ability of lncRNAs.

Figure 3

Specific characteristics of rat lncRNAs. (A) Rat lncRNAs exon number chart. (B) Length distribution of lncRNAs of rat species. (C) lncRNAs orf length statistics. (D) Proportion of lincRNA, antisense lncRNAs, intronic lncRNAs and sense lncRNAs identified in rat.

Figure 4

Expression models of rat lncRNAs. (A) The cumulative curve of the average phastCons score of lncRNAs-exon (blue), lncRNAs-intron (purple), coding genes (red) and genomic background (green). Genomic background was randomly selected from the whole genome. (Graph legend: L01-L03 FPKM: lncRNAs in the normal group, L04-L06 FPKM: lncRNAs in the cold stress group.) (B) The box plot of expression levels of rat lncRNA in the cold stress group and the normal group. The y-axis represents the average log10 (FPKM) value of each set of replicates. (C) Comparison of FPKM density distribution between two groups. (D) Hierarchical clustering of the diferentially expressed lncRNAs.

Figure 5

Differentially expressed lncRNAs. (A) The Volcano plot of DElncRNAs between the cold stress group and the normal group. Red color represents significant up-regulated and green color represents significant down-regulated. (B) The correlation plot, three repeated samples, the points are concentrated in the diagonal, indicating a high correlation coefficient. Correlation was evaluated by Pearson’s correlation coefficient of total lncRNAs expression levels.

Figure 6

GO classification of target genes of DElncRNAs in rats. The abscissa was the GO classification, the left side of the ordinate was the percentage of genes, and the right side is the number of genes.

Figure 7

KEGG classification of target genes of DElncRNAs in rats. The abscissa is the annotated genes, the ordinate was pathway term.

Figure 8

Differential lncRNA interacts with target genes. lncRNAs and coding transcripts were presented as squares and circles, respectively.

Figure 9

The co-expression network of DElncRNAs and DEmRNAs. lncRNAs and mRNAs were presented as diamonds and squares, respectively. And green and blue indicated down-regulate, red indicated up-regulate.

Figure 10

GO classification of co-expressing mRNAs of lncRNA MSTRG.80946.2.

Figure 11

KEGG enrichment analysis of co-expressing mRNAs of lncRNA MSTRG.80946.2.

Figure 12

qRT-PCR validation of high throughput sequencing. Validation of 10 selected DElncRNAs. T-test p-values  < 0.05 are considered to be significantly different, “*” represents a p-value  < 0.05 and “**” represents a p-value  < 0.01.

Figure 13

The agarose gel electrophoresis map of RACE product. M: DL2000 DNA Marker. 5′RAEC: 583 bp (include link sequence); 3′RACE: 383 bp (include link sequence).

Figure 14

The subcellular positioning of MSTRG.80946.2 in BRL cell.

Figure 15

qRT-PCR verification of relative expression after ASO-MSTRG.80946.2 silencing. (A) The relative expression of MSTRG.80946.2. (B) The relative expression of ACP1. (C) The relative expression of Tsn. (D) The relative expression of TSPY1. *P < 0.05, **P < 0.01.