Long noncoding RNAs are involved in multiple immunological pathways in response to vaccination

Significance Long noncoding RNAs (lncRNAs) are known to be involved in several immunological processes. In spite of their general relevance to human immunity, to date there are no reports on the importance of lncRNAs in vaccine responses. Here we apply a “systems vaccinology” framework to study the role of lncRNAs in vaccine-mediated immunity. We applied meta-analytical approaches using public microarray data from over 2,000 blood transcriptome samples of vaccinees and an RNA-sequencing (RNA-seq) dataset from vaccinated children to tackle this question. Our results indicate that lncRNAs are important players in several immunological processes elicited by vaccination.


Quality control and preprocessing of microarray samples
Raw data from microarray datasets were submitted to arrayQualityMetrics package for aberrant sample detection (1). Samples flagged as outliers by two of the three tests performed by the package were removed. The resulting set of samples were iteratively resubmitted to arrayQualityMetrics using the same criteria until no samples were removed. Microarray datasets were then log transformed and quantile normalized. Samples from illumina platforms were normalized using functions provided by limma package (2), whereas samples from Affymetrix platforms were normalized using Robust Multi-Array Average (RMA) from affy (3). Probes targeting genes in common were summarized by selecting the ones showing the highest mean of expression between all samples. Genes whose mean of expression were in the first decile were removed.

Microarray platform reannotation pipeline
We developed a Snakemake pipeline to update probe annotations from microarray platforms. Probe sequences were aligned against the HG38 assembly of the human genome using BLAT (4) with thresholds of alignment score > 90% and identity score > 90%. Probes with multiple alignments were discarded and remaining probes were reassigned to gene feature annotations provided by Gencode v24 using Bedtools (5). Probes aligning with more than one feature were also discarded.

Correlation with immune parameters
Antibody titers, assessed either with Hemagglutination-inhibition assay (HAI) or Microneutralization assay (MN), were retrieved from microarray databases or supplemental materials from articles. Antibody titer fold increases from each subject between day 28, 63 or 70 after vaccination and baseline values were log transformed.
Pairwise log2 fold change values for each day after vaccination (compared to baseline) were correlated (Pearson's correlation) with corresponding antibody titer fold increases.

Correlation between coding and neighboring non-coding transcripts
A catalog of lncRNA-mRNA genomic neighboring pairs was assembled using GenomicRanges package (6). Pairwise log2 fold-change values for each day after vaccination compared to baseline of neighboring pairs were correlated (Pearson's correlation) to infer potential cis-regulatory regions.

Building of a consensus network
Samples collected between days 0 and 7 after IV were submitted to CEMiTool using default parameters (Pearson's correlation coefficient and unbiased selection of genes by a variance-based filter) (7). Genes from each module were connected between themselves to create fully connected subnetworks for each cohort. We then computed the Pearson's correlation coefficient between each pair of genes (edge) in each cohort. The resulting networks were then pruned by selecting edges that: 1) Were represented in at least two cohorts; 2) Had a mean correlation coefficient of greater than 0.6 if edge was represented in two cohorts; 3) Had a mean correlation coefficient of greater than 0.55 if edge was represented in three cohorts; 4) Had a mean correlation coefficient of greater than 0.5 if edge was represented in four cohorts; 5) Had a mean correlation coefficient greater of than 0.45 if edge was represented in five cohorts (only IV); 6) Were represented in 6 or more cohorts for IV and 5 or more cohorts for YF17D. Communities were inferred from resulting consensus network with a spin glass method implemented in igraph R package (https://igraph.org/r/).

Enrichment analysis and Over-Representation Analysis.
Gene Set Enrichment Analysis (GSEA) was performed in pre-ranked lists of individual log2 fold-change values using fgsea package (Bioconductor version: Release 3.8). Over-Representation Analyses were performed using clusterProfiler R package (8). P-values were adjusted, and pathways with False Discovery Rates < 0.05 were selected. Blood Transcriptional Modules (BTMs) defined in (9) were used in both analyses.
Fragments Per Kilobase Million (FPKM) expression data of genes defined in Gencode v24 for cell types matching queries were log transformed.

IV and YF17D vaccination consensus networks
We compared the YF17D and IV networks and found 86 connections between lncRNAs and protein-coding genes that were shared to both vaccines. These connections included FAM30A and the host gene of miR-22 among others. Most connections between lncRNAs and mRNAs were unique to IV network (2,177) or YF17D network (233), suggesting that different lncRNA-mediated regulation may be involved.

GHRLOS
This lncRNA is antisense to ghrelin/obestatin prepropeptide gene (GHRL). The expression of GHRLOS at day 1 post-vaccination is inversely correlated with antibody responses (at day 30) in three IV cohorts. The expression of the sense gene GHRL at day 1 is also inversely correlated with antibody responses (at day 30) in three IV cohorts.
The sense gene Ghrelin plays a critical role in energy homeostasis and it is known as the "hunger hormone". Vishwa et al. (12) showed that ghrelin is expressed in human T cells and it may inhibit proinflammatory cytokines (e.g. TNF, IL-6, IL1B, among others). We suggest that regulation of GHRLOS/GHRL may contribute with the antibody responses induced by IV.

TRIM52-AS1
The expression of lncRNA TRIM52-AS1 at day 3 post-IV is inversely correlated with antibody responses (at day 30) in five IV cohorts (p = 0.00034). Although the TRIM52-AS1 is antisense to TRIM52, the expression of the gene TRIM52 is not correlated with antibody response.
One study shows that down-regulation of TRIM52-AS1 is associated with cell proliferation (13). The Expression Atlas portal (14) shows that TRIM52-AS1 is also down-regulated in Th1 and Th2 polarised T cells (see below). This suggests that the lncRNA may be involved with CD4 T helper differentiation and proliferation.

C1orf132 (MIR29C Host Gene)
The current annotation of lncRNA C1orf132 is MIR29B2CHG which is a host gene of two miRNAs: MIR29C and MIR29B2. The expression of lncRNA C1orf132 at day 7 post-IV is inversely correlated with antibody responses (at day 30) in six IV cohorts (p = 0.0017).
Several studies have shown that the miR-29 family may play a key role in adaptive immunity (15). Chun-Mei et al have shown that miR-29c may suppress the Tumor necrosis factor alpha-induced protein 3 (TNFAIP3), which is a key regulator in inflammation and immunity (16). Since TNFAIP3 is required for the normal differentiation of the marginal zone B and B1 cell subsets in mice (17), we hypothesize that lower levels of miR-29c may be associated with higher antibody-secreting cells by increasing the TNFAIP3 protein levels.

Insights from Influenza vaccination consensus network (Figure 4)
The different cell populations activated upon vaccination may impact the functional status of subsets. For example, the community containing genes associated with antibody-secreting cells (ASCs) is CM7 (Figure 4a (18). In addition to IL1B and NLRP3, CM3 contains a lncRNA (IL10RB-DT) that is antisense to the IL10RB gene, which is essential for the IL-10 receptor complex. Two small nucleolar RNA host genes in CM3, SNHG15 and SNHG17, were previously found upregulated in tumor cells and were associated with cell proliferation (19,20). Their role in CM3 can be related to the proliferation of immune cell types.
Furthermore, monocytes have also been associated with B cell responses and antibody production. In 2014, a research group described that Dengue virus-infected monocytes acquired the CD14+ CD16+ phenotype and were able to elicit the differentiation of B cells into antibody-secreting cells (ASCs) in vitro through the secretion of BAFF (TNFSF13B), APRIL (TNFSF13) and TACI (TNFRSF13B) (21). The ligand BAFF was found in CM6 (which is associated with a type I Interferon response), while its 3 possible receptors TACI, BCMA (TNFRSF17) and BAFFR (TNFRSF13C) were found in CM7 (ASC-related genes). One interesting lncRNA in CM1 is the miR-223 host gene. Several papers have shown that miR-223 is highly expressed in activated monocytes/macrophages and that may regulate NLRP3 inflammasome and IL-1B production (22)(23)(24). Additionally, NLRP12 also belongs to the CM1 module and its expression is negatively regulated by Blimp-1 (PRDM1), which is needed to complete B cell differentiation into ASCs (25). Moreover, Blimp-1 is not only expressed during the final steps of the B cell differentiation into ASCs, but also regulates the differentiation of different T cell subsets (enriched in the CM5) (26,27).