Transcriptomic meta-analysis reveals unannotated long non-coding RNAs related to the immune response in sheep

Abstract

Long non-coding RNAs (lncRNAs) are involved in several biological processes, including the immune system response to pathogens and vaccines. The annotation and functional characterization of lncRNAs is more advanced in humans than in livestock species. Here, we take advantage of the increasing number of high-throughput functional experiments deposited in public databases in order to uniformly analyse, profile unannotated lncRNAs and integrate 422 ovine RNA-seq samples from the ovine immune system. We identified 12302 unannotated lncRNA genes with support from independent CAGE-seq and histone modification ChIP-seq assays. Unannotated lncRNAs showed low expression levels and sequence conservation across other mammal species. There were differences in expression levels depending on the genomic location-based lncRNA classification. Differential expression analyses between unstimulated and samples stimulated with pathogen infection or vaccination resulted in hundreds of lncRNAs with changed expression. Gene co-expression analyses revealed immune gene-enriched clusters associated with immune system activation and related to interferon signalling, antiviral response or endoplasmic reticulum stress. Besides, differential co-expression networks were constructed in order to find condition-specific relationships between coding genes and lncRNAs. Overall, using a diverse set of immune system samples and bioinformatic approaches we identify several ovine lncRNAs associated with the response to an external stimulus. These findings help in the improvement of the ovine lncRNA catalogue and provide sheep-specific evidence for the implication in the general immune response for several lncRNAs.

Keywords: RNA-seq; genomics; immune system; lncRNAs; sheep; transcriptomics.

FIGURE 1

Bioinformatic workflow of the study. The workflow followed in this study can be divided into three sections. (A) First, sequencing data retrieval, preprocessing and mapping to the sheep genome. (B) Second, identification of unannotated lncRNA transcripts and evidence of expression. (C) Third, functional analyses between unstimulated samples and samples with an immune stimulation.

FIGURE 2

Characteristics of the dataset and the identification of lncRNAs. Exploratory analysis of all the samples included in the study using dimensionality reduction methods: (A) Principal Component Analysis (PCA) grouped by main tissue, (B) t-SNE plot with samples colored by tissue. (C) Transcript length distribution of PCGs and lncRNAs. (D) Exon length distribution of PCGs and lncRNAs. (E) Expression levels of PCGs and lncRNAs in blood cell samples and tissue samples. (F) Classification of lncRNAs into classes by genomic location. (G) Number of detected unannotated lncRNAs against sequencing depth.

FIGURE 3

Support for transcription of annotated genes and novel lncRNAs. Fractions of expressed genes with detected TSSs or active gene histone modifications. TSSs were obtained from CAGE-seq peaks from five immune tissues and histone modifications were obtained from ChIP-seq peaks (H3K4me3 and H3K27ac) from alveolar macrophages. PCG: Protein coding gene, Ens_lnc: Ensembl lncRNA, Novel_lnc: Novel lncRNA.

FIGURE 4

Differential expression results between stimulated samples and unstimulated samples. Analyses were performed in blood cell samples (A) and lymph node samples (B). For each comparison, a volcano plot using shrunken fold changes and a dot plot with the results of gene ontology enrichment analysis (GO biological processes) are shown.

FIGURE 5

Co-expression analysis and differential co-expression network results in blood cell samples. (A) Correlations of gene co-expression modules with all stimulations and with each individual stimulation. Modules enriched in immune genes are highlighted in red. Number of genes in each module is depicted as a bar plot. (B) The full differential co-expression network. Node size is proportional to connectivity and differential associations are coloured by gain or loss of correlation strength. The edges of differentially expressed genes are coloured by fold change. (C) Sub-network with the differentially associated genes in module ME16. (D) Sub-network with the differentially associated genes in module ME19.

FIGURE 6

Co-expression analysis and differential co-expression network results in lymph node tissue samples. (A) Correlations of gene co-expression modules with all stimulations and with each individual stimulation. Modules enriched in immune genes are highlighted in red. Number of genes in each module is depicted as a bar plot. (B) The full differential co-expression network. Node size is proportional to connectivity and differential associations are coloured by gain or loss of correlation strength. The edges of differentially expressed genes are coloured by fold change. (C) The genes differentially co-expressed with CREB3 transcription factor. (D) Individual examples of statistically significant differential associations between CREB3 and four genes.