We selected 929 RNA-seq sequencing runs belonging to 15 BioProjects from NCBI Sequence Read Archive (SRA), which were merged into 422 samples, by the following inclusion criteria: Samples from an immune system tissue (blood, immune cells or lymphoid organs), at least five samples from a single BioProject, pair-end sequenced using an Illumina platform and genome mapping rate above 60%. Sample metadata such as tissue type, age, breed, sex, library type or experimental treatment was collected from NCBI databases and published articles. Due to metadata ambiguity, the strandness of the samples was assessed with Kallisto (Bray et al., 2016) prior to pipeline execution.

Most samples originated from functional experiments that studied the immune response to vaccines or vaccine components (PRJEB26387, PRJNA454435, PRJNA559411), helminth infections (PRJNA291172, PRJNA433706, PRJNA268183, PRJEB33476, PRJEB45790, PRJEB44063), bacterial infection (PRJEB15872) and pro-inflammatory gene upregulation (PRJNA631066). Other transcriptomic studies were not related to the immune response but were used to improve the novel lncRNA identification and as unstimulated controls (PRJNA528905, PRJNA485657, PRJNA362606). Besides, we added samples from the sheep expression atlas (PRJEB19199), including samples from bone marrow derived macrophages stimulated with lipopolysaccharide (LPS). All samples were dichotomized into two treatment-groups: samples from immune-stimulated animals and unstimulated or control samples.

We downloaded and analysed the 422 RNA-seq samples with a uniform workflow using custom Snakemake v.6.15.1 (Köster et al., 2021) pipelines (Figure 1). 375 reverse stranded samples were used for transcriptome construction and novel lncRNA identification, while all samples were used for quantification based on the new transcriptome. Sequencing runs were downloaded from NCBI SRA with the SRA Toolkit and were merged into samples by their experiment ID. Adapter trimming and quality filtering was performed with cutadapt v.3.5 (Martin, 2011) using Illumina universal adapters and a phred threshold of 30. Reads were aligned to the sheep reference genome (Oar_rambouillet_v1.0) (Salavati et al., 2020) with STAR v.2.7.3a (Dobin et al., 2013) guided by the Ensembl (v102) annotation. StringTie2 v.2.0 (Kovaka et al., 2019) transcriptome assembler was used to reconstruct the transcriptome of each individual sample guided by the Ensembl (v102) annotation and with the--rf option. Then StringTie2 was applied again with the--merge option using all the transcriptomes in order to obtain a non-redundant transcriptome that is comparable between samples.

Quantification of gene expression was performed at transcript level with Kallisto v.0.48 (Bray et al., 2016) pseudoaligning the trimmed reads of all samples to the newly generated transcriptome, generated with GffRead v.0.11.7 (Pertea and Pertea, 2020). The--rf-stranded option was used with the 375 stranded samples.

Potential novel lncRNAs were defined as unannotated transcripts that were located either in an intergenic region, in an intron of a known gene or in the antisense strand of a known gene. GffCompare v.0.11.2 (Pertea and Pertea, 2020) was used to compare the newly assembled transcriptome with the reference annotation and extract these transcripts. Single-exon transcripts longer than 500 nucleotides and shorter than 10kb, and multiexonic transcripts longer than 200 nucleotides and shorter than 50 kb were kept. The assessment of the coding potential of the candidate transcripts was done with three different tools. The coding potential prediction module of FEELnc (Wucher et al., 2017), based on a Random Forest classifier, was trained with sequences of bovine coding genes and lncRNAs from NONCODE database (Zhao et al., 2021). Coding-Potential Assessment Tool 3.0.2 (CPAT) (Wang et al., 2013) is a logistic regression-based tool that we trained and selected the classification threshold following authors’ instructions using the same bovine coding and non-coding sequences. HMMER 3.3.2 (Eddy, 2011) was used to detect Pfam protein domains in our potential lncRNAs, which were translated into the three possible frames. Transcripts classified as non-coding by FEELnc and CPAT and without protein domains detected by HMMER were kept. Transcripts classified by CuffCompare as a novel isoform of a known gene were also kept, as transcripts that had passed the coding potential tests could be legit non-coding isoforms. The selected transcripts were defined as the final set of novel lncRNAs.

Novel lncRNA transcripts were classified with a custom Python script (see Data Availability section) based on their position relative to their closest gene. Transcription start sites (TSSs) were defined as the start or stop nucleotides, depending on strandness. Seven classes were defined: 1) antisense, for those transcripts overlapping a gene in the opposite strand; 2) sense intronic or antisense intronic, for transcripts fully contained within an intron; 3) intergenic, for lncRNAs at least 5 kb away from any known gene; 4) divergent, with TSSs within 5 kb and in the opposite strand; 5) convergent, with transcription stops within 5 kb and in the opposite strand; 6) sense upstream, located less than 5 kb upstream of a gene and in the same strand; and 7) sense downstream, located less than 5 kb downstream of a gene and in the same strand.

To compare the novel lncRNAs with the recently upgraded ovine NCBI RefSeq annotation (release 104), which is based on the ARS-UI_Ramb_v2.0 new reference genome (Davenport et al., 2022), transcript coordinates were remapped with the NCBI Genome Remapping Service (https://www.ncbi.nlm.nih.gov/genome/tools/remap). They were compared with the NCBI lncRNAs using GffCompare (Pertea and Pertea, 2020). Transcripts models with codes “ = “, “j”, “c”, “k”, “o”,"m" or “n" were considered as overlapping, transcripts with codes “c” or “k” were considered compatible isoforms and transcripts with code “ = ” were considered exact matches.

We downloaded the mapped BAM files of CAGE-seq experiments of five immune tissues from a multi-tissue project of sheep TSSs (tonsil, alveolar macrophages, spleen, mesenteric lymph node and prescapular lymph node) (Salavati et al., 2020) and analysed them using the same pipeline as the authors, with some modifications. In short, downloaded BAM files were converted to bigwig format with bedtools v.2.30.0 (Quinlan and Hall, 2010) and BedGraphToBigWig from UCSC tools (Kent et al., 2010). The R package CAGEfightR v.1.12.0 (Thodberg et al., 2019) was used for normalization and clustering of CAGE tags. CAGE tags <10 read counts were removed and all the remaining tags from any of the tissues were kept, to include tissue-specific TSSs. CAGEfightR was also used to identify bidirectional clusters. In order to get the genes supported by CAGE-predicted TSSs we used the BedTools python implementation pybedtools v.0.8.1 (Dale et al., 2011) to search for TSSs from the assembled transcriptome within 0.5 kb from them, accounting for strandness.

Sheep ChIP-seq sequencing files from alveolar macrophages (Massa et al., 2021) were downloaded from the NCBI Sequence Read Archive (SRA) and were analysed in an uniform way. Reads were aligned to the sheep genome (Oar_rambouillet_v1.0) with Bowtie2 v.2.3.5.1 (Langmead and Salzberg, 2012). SAM files were converted to BAM format with samtools v.1.7 (Danecek et al., 2021), and were sorted, filtered for quality and removed duplicate reads with sambamba v.0.6.6 (Tarasov et al., 2015). MACS2 v.2.2.6 (Zhang et al., 2008) was used to call narrow peaks for histone modifications with a FDR cut-off of 0.05 and consensus peaks from the pairs of animals were obtained with bedtools v.2.30.0 (Quinlan and Hall, 2010). In order to get the genes supported by ChIP-seq peaks we used pybedtools v.0.8.1 to search for TSSs from the assembled transcriptome within 0.5 kb from them.

Sequence level conservation was performed with standalone BLASTn (BLAST v.2.9.0) (Camacho et al., 2009) by aligning the sheep lncRNA transcripts against the lncRNAs annotated in Ensembl Release 106 from five species: goat, cattle, pig, mouse and human. Because of the known low sequence conservation expected in lncRNAs, results were filtered by identity >50, query coverage >50, E-value > 1e-05 and it was required that the length differences between each pair of sequences was less than 50%. Visualization of the genomic context of conserved lncRNAs was performed with pyGenomeTracks 3.7 (Lopez-Delisle et al., 2021). The tracks for CAGE-seq data were constructed by merging all BAM alignment files with samtools (Danecek et al., 2021) into a single file and then was converted to bigwig format as previously. The tracks for histone modification ChIP-seq data were the consensus peaks obtained from MACS2.

Kallisto abundance estimates were imported to R and summarized to gene level with IsoformSwitchAnalyzeR (Vitting-Seerup et al., 2019) in order to set confident gene identifiers for ambiguous transcripts. Counts of annotated genes and novel lncRNA genes were kept for further analysis, discarding potential novel unannotated coding genes. For gene expression data exploration, we normalized the estimated gene counts with the variance stabilizing transformation from DESeq2 (Love et al., 2014) and filtered out genes with less than 0.5 transcripts per million (TPM). The first two components of the principal component analysis (PCA) and the two first dimensions of the t-Distributed Stochastic Neighbor Embedding (t-SNE) were used for visualization. LncRNAs were tagged as expressed if they could be detected above 0.1 TPM or one TPM in at least 20% of the samples in a tissue group. Two-sided Mann-Whitney U tests were performed to compare expression means between classes.

Differential gene expression (DGE) was performed with DESeq2 (Love et al., 2014) using the estimated counts of annotated genes and lncRNA genes expressed in at least half of each sample groups and exported from IsoformSwitchAnalyzeR (Vitting-Seerup et al., 2019). Differential expression was tested separately in a blood and cell sample dataset on one side, and in a lymph node-only dataset on the other, because there were not stimulated samples from other lymphoid organs and that would unbalance the dataset. The Wald test was applied between unstimulated samples and stimulated samples using the effect of the interaction of tissue type and BioProject IDs as covariates for the lineal regression model, as those were the main drivers of the groupings seen in the exploratory analysis. Log₂ fold change (log₂FC) values from lowly expressed and highly variable genes were shrunken using the apeglm method (Zhu et al., 2019). Genes with an FDR-adjusted p-value lower than 0.05 and an absolute log₂FC higher than 0.32, which corresponds to a 20% expression change, were kept. The relatively low log₂FC filter was chosen because the large number of samples and the heterogeneity of the dataset produced differentially expressed genes with modest effect sizes and robust p-values.

Gene set enrichment analysis of differentially expressed genes was done with gProfiler R package (Raudvere et al., 2019). The statistical domain scope used was the list of all expressed genes for each tissue, in order to reduce the tissue type specific expression bias. Benjamini–Hochberg FDR correction was applied to the p-values and the threshold was set to 0.05.

Co-expression analyses were performed in the blood and immune cell dataset and in the lymph node-only dataset separately. Genes expressed in less than half of the samples were removed and strong outlier samples were removed in order to get a better fit to a scale-free topology. We tested two network construction pipelines: 1) The pipeline proposed by the authors of Gene Whole co-Expression Network Analysis (GWENA) (Lemoine et al., 2021), which consists of applying the variance stabilizing transformation (VST) from DESeq2 (Love et al., 2014) and using spearman correlations, and 2) counts adjusted with trimmed mean of M-values (TMM) factors followed by asinh transformation, Pearson correlations and network transformation by context likelihood of relatedness (CLR) (Johnson and Krishnan, 2022). Before creating the correlation matrices, normalised gene expression was corrected for covariates with limma’s removeBatchEffect function (Ritchie et al., 2015) to account for the effect of the interaction of tissue type and BioProject ID, as those were the main drivers of the groupings seen in the exploratory analysis. The 30% less variable genes were removed for network construction. Co-expression networks were constructed with GWENA (Lemoine et al., 2021) R package, which implements the Weighted Correlation Network Analysis (WGCNA) (Langfelder and Horvath, 2008) R package.

Modules of co-expressed genes were detected with the threshold power and clustering threshold calculated by GWENA and a minimum module size of 30. Modules were merged if their eigengene, the first principal component of the module, correlation was higher than 0.9. Modules were associated with overall immune stimulation or specific stimulation types by correlating their eigengene to those variables. To calculate the correlation p-value threshold, we generated 1,000 random gene modules ranging from 30 to 1,000 genes, correlated their eigengenes with the treatment variable and calculated the false positive rate (FPR). The p-value threshold with the FPR lower than 0.05 was 1e-02, but 1e-03 was chosen for more robustness. The genes in each module were tested for Gene Ontology (GO) term enrichment with gProfiler (Raudvere et al., 2019) R implementation, setting the statistical domain scope to all the genes in the co-expression network and a FDR-adjusted p-value threshold of 0.05.

The differential co-expression analysis was carried out by calculating the spearman correlations between all genes used in the co-expression network analysis separately in the unstimulated and the stimulated samples. The z-score method implemented in the dcanr v.1.12.0 R package (Bhuva et al., 2019) was used for testing the statistical differences between z-transformed correlation coefficients in both conditions. p values were adjusted for multiple hypothesis testing in order to select differentially correlated gene pairs. Differential co-expression networks (DCN) were visualized in Cytoscape v.3.8.2 (Shannon et al., 2003) by integrating the differential co-expression results, co-expression modules and differential expression results. For visualization, genes without gene names in the Ensembl annotation were named after their human orthologue according to Ensembl Compara.

We collected and analysed 422 publicly available RNA-seq samples of tissues related to the immune system (Supplementary Table S1) using a uniform pipeline (Figure 1). In terms of immune response induction, 49.1% of the samples had been stimulated in some way. Blood samples, as whole blood or PBMCs, represented the 64.5% of the dataset, organs and lymph nodes the 30.3% and immune cell subsets the 5.2%. The mean age of the animals was of 1.32 years and 60.1% of the samples came from male sheep. There are 12 different breeds in the dataset, with three of them being crossbreed. Library selection is an important factor for lncRNA profiling because there are transcripts that are not polyadenylated. Around half of the samples were polyA-selected and half of the samples sequenced total RNA. Besides, samples had an average of 45 million reads, summing around 19 billion reads in total (Supplementary Figure S1). Unique genome mapping rate with STAR was of 84.6% on average and pseudoalignment rate to the new transcriptome with Kallisto was of 84.7% on average (Supplementary Figure S1). The assembled and merged transcriptome annotation contained 308750 transcripts, of which 41638 were from annotated transcripts and 36067 were from novel lncRNA transcripts identified by our pipeline. These transcripts correspond to 63364 genes, including 25472 annotated genes and 21223 novel genes with at least one lncRNA isoform.

All samples were clustered based on gene expression to assess the coherence of the data. Both clustering methods used clustered together the samples based on tissue, although intra-tissue groupings were influenced by the source project (Figures 2A,B). This could be expected as each study was performed in different conditions, with different breeds, ages, sex and protocols. Immune stimulation status did not affect much the clustering probably for the same reasons and because of the strong influence of tissue type.

We identified 21223 unannotated lncRNA genes from the sheep immune system samples that were assembled, and another 1724 annotated genes had novel non-coding isoforms classified as lncRNAs by our pipeline. Most of the novel genes with transcripts fulfilling the requisites to be classified as novel lncRNAs had all of their isoforms classified as such (17605). Some of the newly assembled gene models were coding genes missing from the Ensembl annotation that had non-coding isoforms, because they had novel transcripts with coding potential as well as lncRNA transcripts. Those unannotated genes and the 1724 annotated genes with novel non-coding isoforms were not considered as lncRNA genes for the gene-level expression analyses, even if individual transcripts could not be discarded as bona fide non-coding isoforms. We applied the same coding potential assessment methods used for novel transcripts to the annotated lncRNAs and discovered that there were transcripts potentially coding for a protein. This bias should be taken into account when comparing between the features of annotated and unannotated lncRNAs.

Regarding the characteristics of the novel lncRNAs, novel transcripts were shorter than the 2,229 lncRNAs annotated in Ensembl (Figure 2C) and a great proportion of them had two exons, in contrast to the Ensembl lncRNAs, which are monoexonic or have more than five exons (Figure 2D). We classified the novel genes based on position relative to known genes (Figure 2F, Supplementary Table S1). Intergenic lncRNAs (lincRNAs) were the most prevalent with 37% of the transcripts, followed by intronic antisense (22%) and antisense (18%) transcripts. Among the transcripts adjacent to annotated genes, the class of divergent lncRNAs was predominant (10% of all novel genes). The TSS of this kind of lncRNAs are very close to another gene’s TSS, which indicates that they probably arise from a single bidirectional promoter and may have implications in terms of gene expression regulation.

We explored the sequence-level evolutionary conservation of lncRNAs with other mammal species. Most lncRNAs are known to be poorly conserved in terms of sequence, but by detecting mammalian orthologues we provide further strength to the methods by which all unannotated lncRNAs have been identified. This analysis found a small number of conserved lncRNAs (Supplementary Figure S2; Supplementary Table S3). The biggest fractions of lncRNAs with conserved orthologues were found when comparing with goat and cattle lncRNA catalogues, with 11.9% and 7.3% of transcripts with significant hits, respectively. Comparing with the human and mouse catalogues, we found much less conserved lncRNAs. Interestingly, around 3% of novel lncRNAs, corresponding to 746 unique transcripts, matched with 392 unique human lncRNAs. Among these conserved lncRNAs, widely characterized lncRNAs such as MALAT1, NEAT1, XIST, PACERR or FIRRE were successfully detected in sheep (Supplementary Figure S3). Other conserved lncRNAs were those located in Hox gene loci, such as HOTAIR, HOXA10-AS, HOXA-AS2 or HAGLR. Divergent lncRNAs were also among the conserved ones, like FMNL1-DT, TOB1-AS1, EMSY-DT, RIPK2-DT, ATP8A1-DT or MAPK6-DT. Despite not showing enough sequence similarity, we found some sheep transcripts located in the same divergent promoter as their human counterparts, for instance the putative orthologues of HEATR6-DT or NIPBL-DT.

Because of the recent improvement of the ovine NCBI reference genome and annotation (Davenport et al., 2022), the NCBI RefSeq lncRNA annotation was compared with the novel lncRNAs. After remapping to the new genome, out of the 4442 NCBI lncRNA genes, 1961 (44%) overlapped with an unannotated lncRNA. Exact matches of intron chains occurred in 571 transcripts, 238 transcripts were intron-compatible but differed in exon number and 3,679 where multi-exonic transcripts with at least one intron match. The overlap between Ensembl and NCBI lncRNAs was virtually inexistent. Thus, we detected around half of the annotated NCBI lncRNA genes using only immune-related tissues, even if most of the transcript models diverged in terms of splice-junctions.

Expression levels of the novel lncRNAs detected in this study were lower than both protein coding genes and other annotated lncRNAs in the two main tissue categories (Figure 2E). In fact, after applying a minimum expression threshold in each tissue, expressed in at least 20% of the samples with one TPM, we were left with 2,267 expressed novel lncRNAs. Besides, we also detected 482 annotated lncRNAs above the expression threshold. Interestingly, 70% of the lncRNAs annotated by Ensembl were expressed in all three main tissue categories, while only 15% of novel lncRNAs were expressed in the three tissues (Supplementary Figure S4). Setting a less stringent mean expression threshold of 0.1 TPM results in 10045 expressed novel lncRNAs, 28% of them in all three tissues. Most of the novel lncRNAs (87%) and annotated lncRNAs (93%) could be detected in the set of lymphoid organs. Overall the overlap was greater between the blood samples and “immune cell” samples for both lncRNA genes and protein coding genes, as blood contains most of those cells (Supplementary Figure S4).

The amount of detected lncRNAs in each sample significantly correlated with sequencing depth for both unannotated lncRNAs (Pearson r = 0.75) and annotated lncRNAs (Pearson r = 0.85). Expression of PCGs was also correlated (Pearson r = 0.58) with sequencing depth but the saturation curve showed a flatter slope, meaning that it saturated earlier than lncRNAs (Supplementary Figure S5). The amount of lncRNAs expressed above 0.1 or one TPM got saturated above around 50 million reads, while the overall amount of expressed lncRNAs at any level did not saturate even at the highest sequencing depths in the dataset (above 100 million reads) (Figure 2G).

Divergent lncRNAs showed greater expression levels than other lncRNAs classes such as intergenic lncRNAs (Mann-Whitney U test p-value 2.9e-10) or antisense lncRNAs (Mann-Whitney U test p-value 1.3e-03), and only showed significantly lower levels than convergent lncRNAs (Mann-Whitney U test p-value 2.4e-03) (Supplementary Figure S6). Intronic antisense lncRNAs showed consistently lower expression than the rest of novel lncRNAs classes, in contrast with convergent lncRNAs, which were significantly more expressed than all other classes.

Independent datasets of CAGE-seq and histone modification ChIP-seq were used in order to provide evidence of lncRNA transcription at RNA and DNA level. The CAGE-seq dataset contained samples from various lymphoid organs and alveolar macrophages, so it was used to provide support of expression in two sample subsets, blood and other immune cells, and lymphoid tissues. We obtained over two million significant CAGE peaks and around 30 thousand bidirectional CAGE peaks present in any of the five tissues.

In both sample subsets PCGs were more strongly associated with CAGE peaks than lncRNA genes, but reducing the analysis to the genes expressed above one TPM instead of 0.1 TPM increased the support in all gene types (Figure 3). This increase in support specially happened in lncRNAs. 64% and 50% of the TSSs of novel lncRNAs expressed above one TPM in the blood subset and the lymphoid subset, respectively, were located within500 bp of a CAGE peak. LncRNAs annotated by Ensembl reached a support level comparable to that of PCGs, with more than 90% of supported TSSs at one TPM. Bidirectional CAGE tag clusters are usually used to identify active enhancers because it is known that bidirectional transcription of short transcripts, known as enhancer RNAs (eRNAs), is a hallmark of enhancer activation (Andersson et al., 2014). Considering the genes expressed above one TPM or 0.1 TPM, novel lncRNAs were slightly less enriched in bidirectional clusters than PCGs. Around 11% and 8% of novel lncRNAs in the blood and tissue datasets, respectively, were transcribed from bidirectional sites (Figure 3). Some of them could be enhancer associated non-coding transcripts while others are divergent lncRNAs.

As for the ChIP-seq data, we analysed two histone modifications that are relevant for lncRNA transcription from a published dataset: H3K4me3, associated with promoters (Santos-Rosa et al., 2002), and H3K27ac, associated with active enhancers and promoters (Creyghton et al., 2010). The trend of H3K4me3 peaks from alveolar macrophages and CAGE peaks were similar regarding the genes expressed in blood and other immune cells, but the overlap between histone ChIP-seq data with TSSs randomly located in the genome was much lower (Figure 3). PCGs had the highest proportion of these promoter-associated marks followed by annotated lncRNAs and novel lncRNAs. Nevertheless, regarding the H3K27ac modification, the difference between lncRNAs and PCGs was smaller, which reflects the origin of many lncRNAs from enhancer-like regions. The support from this modification was similar in novel lncRNAs and annotated lncRNAs. 20% of the TSSs of novel lncRNAs expressed above 0.1 TPM in the blood subset were associated with H3K27ac. The apparent higher support for annotated non-coding models is probably linked with their misannotation.

Providing additional evidence of the transcription of novel transcripts assembled from short-read RNA-seq reads ensures that the detected genes are reproducible. We selected 12302 assembled gene models as bona fide lncRNA genes, those which were supported by at least one of the following: CAGE tags, histone modification ChIP-seq peaks or expressed above 0.1 TPM in at least 20% of the samples in a tissue group (Supplementary Figure S7). In this set, 47% of the lncRNAs had at least support from CAGE peaks or histone modifications. Around 1,000 lncRNAs were supported by all assays, including both histone modifications. The annotation files with all unannotated lncRNA transcripts, the set of high confidence transcripts and expression values can be found in a public repository (Bilbao-Arribas and Jugo, 2022).

We performed differential expression analysis between unstimulated or control samples and samples stimulated with either vaccines or a pathogen in order to identify common lncRNAs induced during an immune response. In blood samples there were 716 differentially expressed genes, including 75 novel lncRNAs and 22 annotated lncRNAs (Figure 4A; Supplementary Table S4). The large number of samples used in the blood sample dataset 222) and the heterogeneity of the data produced many differentially expressed genes with modest effect sizes but robust p-values (Supplementary Figure S8). The most significant enriched terms among the known genes were biological processes related to the immune response to external pathogens such as response to external stimulus (GO:0009605, FDR = 2.86e-09), response to virus (GO:0009615, FDR = 6.75e-07) or defense response (GO:0006952, FDR = 1.19e-07). In lymph node samples, there were 365 differentially expressed genes, including 46 novel lncRNAs and 13 annotated lncRNAs (Figure 4B; Supplementary Table S4). In this case, among the most significant enriched terms with the highest quantity of genes were general terms such as response to stress (GO:0006950, FDR = 2.51e-04) and response to stimulus (GO:0050896, FDR = 1.21e-03). More specifically, the terms related to T cell activation, like T cell activation (GO:0042110, FDR = 3.07e-03) and regulation of T cell activation (GO:0050863, FDR = 3.37e-03), reflect the critical roles of lymph nodes in adaptive immunity. Besides, there also were highly significant but smaller in size enriched terms related to response to endoplasmic reticulum (ER) stress.

There were 22 differentially expressed genes common to both datasets, among them an annotated lncRNA and an unannotated lncRNA. Some of the common PCGs are directly related with immunity, like IL21, which encodes a well known cytokine with immunoregulatory activity that induces proliferation and differentiation in several immune cell types. Other genes are related to apoptosis and inflammation (MT2, IKBIP, AEN, OSGIN1) and ER regulation (WFS1, SELENOS). Despite relatively similar number of DE genes in both comparisons, there is a big set of highly significant genes with effect sizes smaller than the threshold in the blood samples and many statistically significant but lowly expressed novel lncRNAs did not pass the fold change threshold because they were shrunken (Supplementary Figure S8). These results give support for potential involvement of a fraction of the detected novel and annotated lncRNAs in both the innate and adaptive immune responses, following the guilt-by-association principle.

Gene co-expression networks were constructed providing valuable information about the expression relationships of lncRNAs with PCGs and allowing the inference of their putative functions by guilt-by-association. We tested two different network construction pipelines and selected the one proposed by the authors of GWENA (Lemoine et al., 2021), as it produced networks with better fit to a scale-free topology and most of the genes could be associated to an expression module. Covariate correction for tissue type and source project enabled the construction of unbiased networks (Supplementary Figure S9.10). Filtering of lowly expressed genes, genes with low variability and outlier samples that reduced the fit to a scale-free topology resulted in co-expression networks of 12898 and 13428 genes in blood samples and lymph nodes, respectively. In the blood dataset, genes with similar expression patterns were clustered in 33 modules ranging from 54 to 1832 genes (Figure 5A; Supplementary Table S5), and in the lymph node dataset genes were clustered in 30 modules ranging from 44 to 1909 genes (Figure 6A; Supplementary Table S5). Most modules included novel lncRNAs and annotated lncRNAs, and some of them were even hub genes of their module.

We searched for significant correlations among module eigengenes, the principal component of the genes in the module that depicts its dominant trend, and treatment variables. In the blood sample dataset, 15 modules were correlated (p-val < 1e-03) with the general treatment variable, which accounts for any kind of sample stimulation (Figure 5A). Considering correlations to specific immune stimulations, helminth infection shared many correlated modules with the general treatment variable, which meant that it was one of the main drivers of variability in the dataset. Stimulation with LPS and with inactivated foot-and-mouth disease virus (iFMDV) were correlated with specific gene modules different to those correlated to helminth infection. Other stimulations were also correlated to some modules but because of their small sample size they were not further taken into account.

Gene expression modules were characterised by GO term enrichment (Supplementary Table S5). Two of the stimulation-correlated modules (ME16 and ME19) were highly enriched in biological processes related to the immune response but they were not correlated with helminth infection. ME16 was associated with the sum of all treatments and was specially strongly correlated with LPS stimulation. The most significant enriched biological process GO terms were related to the general immune response, like immune system process (GO:0002376, FDR = 2.27e-07) or immune response (GO:0006955, FDR = 1.15e-05) and to cell migration and locomotion, including the terms positive regulation of locomotion (GO:0040017, FDR = 1.55e-06) and leukocyte migration (GO:0050900, FDR = 1.16e-05). ME19 was also associated with the sum of all treatments and was correlated to iFMDV treatment. It contained a high amount of immune response genes, for instance, from the 155 genes with GO annotations, 35 were related to response to virus (GO:0009615, FDR = 9.59e-25) and 56 to immune system process (GO:0002376, FDR = 1.56e-10). Besides, terms related to type I interferon response and signalling were also abundant.

In the lymph node network, the eigengenes of 11 modules were correlated (p-val < 1e-03) with the combined treatment variable (Figure 6A). The two available immune stimulation conditions, helminth infection and paratuberculosis, were correlated with a few modules, but several other significant modules emerged from the combined treatment variable correlation. The characterisation of gene expression modules by GO term enrichment revealed up to five immune-enriched modules: ME15, ME19, ME24, ME27 and ME28. Among them, the positively correlated modules showed functions involved in the innate immune response and general immune terms. For instance, in module ME15 the terms immune response (GO:0006955, FDR = 2.10e-11) or innate immune response (GO:0045087, FDR = 5.77e-07) are highly significant. In contrast, the negatively correlated modules are enriched in adaptive immune response terms. ME19 is enriched in GO terms related with T cell activation and lymphocyte proliferation while ME27 is enriched in terms related to B cell activation and proliferation. The lncRNAs present in the immune-enriched modules from both co-expression networks were classified as immune response-related lncRNAs.

In addition to the immune-enriched gene modules, another big module stood up (ME3), as it was correlated with both helminth infection and paratuberculosis. Most of the enriched GO terms were related to endoplasmic reticulum (ER) stress and protein post-translational processing, with terms like response to endoplasmic reticulum stress (GO:0034976, FDR = 2.60e-13), Golgi vesicle transport (GO:0048193, FDR = 2.42e-07) or response to unfolded protein (GO:0006986, FDR = 1.72e-06).

Gene level differential co-expression, the gain or loss of correlation between two genes in different biological situations, indicates changes in regulatory relationships between those genes, which are often not evident from DGE results. All gene-pairs used in the co-expression network construction were tested for significant changes in correlation between control and stimulated samples and differential co-expression networks (DCN) were constructed with statistically significant gene-pairs (Supplementary Table S5). The DCN from the blood sample dataset contained 1,589 differential associations (FDR <0.05) among 1,348 genes (Figure 5B) and the DCN from lymph nodes contained 2,137 differential associations (FDR < 1e-03) among 1784 genes (Figure 6B). Both networks included around 60 lncRNAs each. In terms of network topology, networks showed a small amount of nodes (genes) with many edges (differential associations), while the rest of the nodes were more loosely connected to the network. Just around 5% of the nodes had 10 or more edges. Some very interconnected nodes formed clusters according to the gene co-expression modules from the previous analysis, but most of the topology was driven by a few high-degree nodes.

Specific differential associations were observed by individually inspecting each DCN. The blood sample network was centred on two high-degree genes that had more than 100 differential associations each but did not have obvious biological relationship with the immune response: DNAJB4 and GUCY1B1. Among the rest of the 42 high-degree genes, defined as those with more than five differential associations, there were some lncRNAs and several immune genes such as BATF2, IDO1, IFI6, IL18BP, NFKBIZ and various CC chemokines. Focusing on the immune-enriched co-expression modules, many genes from the module ME16, most of them immune-related, formed a very interconnected subnetwork (Figure 5C). Even though the genes from this subnetwork were already correlated, they predominantly showed positive z-scores, which means that the correlations were stronger in the control samples than in stimulated samples, On the contrary, genes from the module ME19 did not form a separate cluster, but they showed negative z-scores, which means that their expression was correlated in the stimulated samples. Interestingly, many genes were up-regulated in the differential expression analysis. In the subnetwork composed by selecting the genes from this module and their differentially co-expressed pairs, there were transcription factor coding genes related to the immune response: BATF2, IRF9 and NFKB2 (Figure 5D). For instance, BATF2, upregulated in stimulated samples, is a transcription factor that controls the differentiation of lineage-specific cells in the immune system and immune-regulatory networks. There were several interferon-stimulated genes such as IFI6, MX1, MX2, ADAR, EIF2AK2, IRF9 or IFIH1, all related to antiviral functions and upregulated in the stimulated samples.

The DCN obtained from lymph node samples did not contain many immune-related genes (Figure 6B). There were 187 high-degree genes, including 18 transcription factors coding genes that were potential drivers of the differential co-expressions, like HMBOX1, CREB3, NFATC4, NFIB or EBF4. NFATC4, for instance, is involved in T-cell activation, stimulating the transcription of IL2 and IL4 cytokine genes. CREB3, among many other functions, plays a role in the response to ER stress by promoting cell survival, a process that was previously found enriched in a co-expression module of which CREB3 was not part of. CREB3 was a high-degree node, differentially associated with 27 other genes (Figure 6C), and it showed mostly positive z-scores, thus, its expression was correlated in the unstimulated samples but those correlations were lost upon stimulation by helminth infection and paratuberculosis. Examples of differentially associated genes include IL18BP, CTSL, MAPK13 and ADM (Figure 6D). ADM, which codes for a known lymphangiogenic factor, was upregulated in stimulated samples and its expression decoupled from that of CREB3 in those samples. This DCN also contained several lncRNAs and 12 of them were high-degree nodes (7 known lncRNAs and 5 novel lncRNAs).

We used differential gene expression analysis, co-expression analysis and differential co-expression network analysis for the functional association of lncRNA genes with the activation of the immune response. Those three approaches resulted in 320 lncRNAs associated in at least one analysis (Supplementary Figure S11). The differential expression between stimulated and unstimulated samples showed the highest number of immune response-associated lncRNAs. Interestingly, the histone modification support in differentially expressed novel lncRNAs was much higher than in the whole set of novel lncRNA genes, 49% against 19%, and the trend was similar in the case of CAGE support. A summary of all transcription evidence and associations in an analysis for each lncRNA is available as a supplementary file (Supplementary Table S6).