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

Vaccine Cohorts and Reannotation of Microarray Platforms.

Blood and peripheral blood mononuclear cell (PBMC) expression data from 17 cohorts of human participants immunized with either yellow fever vaccine (YF-17D) or inactivated influenza vaccine (IV) were analyzed (Fig. 1A and SI Appendix, Table S1). Most of these cohorts included antibody titer data collected before and after vaccination (Fig. 1A). RNA derived from these samples was hybridized onto microarrays, and raw data were deposited in the Gene Expression Omnibus (GEO) and ArrayExpress databases. After data preprocessing, four types of analyses for each dataset were conducted to identify lncRNAs associated with vaccine-induced immune responses (Fig. 1B and Methods).

• Download figure
• Open in new tab
• Download powerpoint

Fig. 1.

Experimental design of meta-analysis of 17 vaccine cohorts. (A) Graphic representation of sampling regimen. Blood samples were collected on days 0, 1, 3, 7, 14, and 28 after vaccination with IV or YF-17D. Numbers in circles represent the number of samples used from each cohort after processing. The Influenza season is shown in the name of the IV cohort (07 = 2007/2008, 08 = 2008/2009, etc.). (B) Summary of computational analyses performed in this work. Blood transcriptome from cohorts was assessed with commercial microarray platforms. After remapping of probes, each cohort was subjected to analysis of differential expression, correlation to antibody titers, correlation to neighboring genes, and coexpression analysis. (C) Characterization of lncRNA classes represented in each microarray platform.

To identify probes representing lncRNAs, the probe sequences from each microarray platform were aligned to the human genome assembly hg38 and were reannotated based on genomic features described in GENCODE version 24. All microarray platforms contained probes targeting different classes of lncRNAs (Fig. 1C and SI Appendix, Fig. S1 A and B). Some probes that were officially annotated by the microarray manufacturer as representing protein-coding genes were, in reality, targeting antisense, intronic, or intergenic lncRNAs (SI Appendix, Fig. S2). Since each microarray platform contains probes targeting specific sets of genes and lncRNAs, we sought to characterize intersections between platforms. A total of 204 and 91 lncRNAs were detected simultaneously in all microarray platforms used in IV and YF-17D vaccine cohorts, respectively. These numbers increased when different combinations of microarray platforms were considered (SI Appendix, Fig. S1 A and B).

lncRNAs Present Coherent Changes of Expression after Vaccination.

Dozens of lncRNAs were identified whose expression was induced or repressed after 1, 3, 7, or 14 d following IV compared with baseline (Fig. 2A and SI Appendix, Fig. S3). Among the protein-coding genes, TNFRSF17, GGH, and CD38 were up-regulated 7 d postvaccination in most cohorts (Fig. 2 B and C), as previously reported (6). Gene set enrichment analysis (GSEA) was performed for each comparison using blood transcription modules (BTMs) (26). For IV, signatures related to monocytes, cell cycle, TLR signaling, antigen presentation, and IFN response were mostly induced on day 1 postvaccination, whereas a strong signature associated with B cells and stimulated CD4+ T cells were detected in the majority of cohorts on day 7 postvaccination (SI Appendix, Fig. S4).

• Download figure
• Open in new tab
• Download powerpoint

Fig. 2.

Transcriptome analysis of cohorts immunized with inactivated IV. Cumulative sum of differentially expressed lncRNAs (A) and protein-coding genes (B) (y-axis) in one or more datasets (x-axis). Transcripts were considered differentially expressed if limma P values were lower than 0.05 in at least three cohorts (for lncRNAs) or four cohorts (for protein-coding genes). (C) Forest plots of representative lncRNAs and TNFRSF17 with reiterated differential expression. Log₂ fold changes with their corresponding 95% confidence interval (x-axis) are plotted for each cohort (y-axis). Red vertical lines and shaded regions represent log₂ fold-change summaries and their 95% confidence intervals, respectively. (D) Heatmap depicting gene expression of human immune cells from ref. 28. Logs (FPKM) of genes are scaled around zero. Columns represent samples; rows represent genes.

We also performed random-effects meta-analyses (27) to identify lncRNAs with coherent differential expression throughout the IV and YF-17D cohorts. Often, lncRNAs showed subtle but significant changes in expression following vaccination. For the IV cohorts, MIAT and MCM3AP-AS1 were down-regulated on day 1 postvaccination (Fig. 2C) (nominal P values = 3.05 × 10⁻⁶ and 1.06 × 10⁻⁵, respectively). LINC01133 and DANCR were up-regulated on day 3 (nominal P value = 0.000917) and on day 7 (nominal P value = 1.05 × 10⁻⁸) postvaccination, respectively (Fig. 2C). As expected, the protein-coding gene TNFRSF17 was up-regulated on day 7 post-IV in all cohorts (nominal P value = 1.787 × 10⁻¹¹) (Fig. 2C). For YF-17D–vaccinated cohorts, differentially expressed lncRNAs included AATBC (nominal P value = 0.0008938 on day 7), MAPKAPK5-AS1 (nominal P value = 0.0005284 on day 7), and LRRC75A-AS1 (nominal P value = 0.000339 on day 3 and 2.524 × 10⁻⁸ on day 7). Increases in the expression of IFN-stimulated genes after vaccination is a mark of YF-17D response (3, 4). Accordingly, we detected sharp and coherent increases in the expression of IFI27 in all YF-17D cohorts on day 7 (nominal P value = 6.393 × 10⁻⁴⁵) (SI Appendix, Fig. S5).

Next we investigated the expression of these lncRNAs in different immune-related cells using RNA-Seq data from Blueprint (28) (see Methods, Fig. 2D, and SI Appendix, Fig. S4). DANCR, which had previously been implicated in the control of histone methylation mediated by EZH2 (29, 30), is specifically expressed in cells of lymphoid origin (Fig. 2D). The lncRNA AATBC, which is associated with cell proliferation and apoptosis (31), was expressed specifically in monocytes (SI Appendix, Fig. S5).

FAM30A lncRNA Levels Correlate with Those from Neighboring Immunoglobulin Heavy Genes.

We also identified lncRNAs whose expression correlated with antibody response in both IV and YF-17D cohorts. The lncRNA with the highest correlation throughout all IV cohorts was FAM30A (KIAA0125), a lncRNA embedded in the immunoglobulin (Ig) heavy (IgH) locus on chromosome 14, in antisense orientation to IgH gene segments (Fig. 3). The expression of FAM30A showed a positive correlation to the levels of antibody titers on day 7 postvaccination in most of the IV cohorts (Fig. 3A). Higher expression levels of FAM30A were seen in B cells compared with other immune cell types (Fig. 3B). Increases in the frequency of antibody-secreting cells (ASCs) could be observed 7 d post-IV (11, 13) and may partly explain the observed positive correlation between FAM30A expression and antibody titers seen in the microarray analysis. To check for evidence of coregulation, we compared the expression levels of FAM30A with those of 4 IgH locus genes located in its genomic vicinity on day 7 post-IV. This approach revealed a strong positive correlation with IgH gene segments (Fig. 3C), suggesting potential related functions and a cis-regulatory role for FAM30A in regulating the expression and functional organization of the gene segments. We describe examples of lncRNAs whose expression is inversely correlated with antibody responses in the SI Appendix.

• Download figure
• Open in new tab
• Download powerpoint

Fig. 3.

FAM30A expression correlates to antibody production and to gene segments within the IgH locus. (A) Correlations between FAM30A and antibody titers with corresponding 95% confidence interval (x-axis) are plotted for each cohort (y-axis). Red vertical line and shaded region represent the correlation summary and its 95% confidence interval, respectively. (B) FAM30A expression is higher in B cells according to ref. 28. (C) Fold-change correlation between FAM30A and gene segments within the IgH locus at day 7 after inactivated IV. Each circle represents Pearson correlation coefficient in different IV cohorts. Meta-analysis P values are represented below cohort identifiers.

Modular Expression Analyses of lncRNAs.

Assessing gene coexpression patterns in blood transcriptome may help in the identification of novel functional connections between lncRNAs and mRNAs after vaccination (32). To detect shared coexpression patterns between IV cohorts, we first ran CEMiTool (33) in each cohort using samples from baseline to day 7 after vaccination. Next we integrated the resulting modules into a consensus network and used a spin glass method to identify 16 cohesive network subgroups (34), hereafter referred as communities (Fig. 4A). These communities were enriched for distinct immune-related cells and processes such as monocytes (CM1), platelet activation (CM4), T cells (CM5), IFN response (CM6), B cells (CM7), and natural killer cells (CM9) (Fig. 4B). GSEAs were then undertaken using the communities as gene sets and mean log₂ baseline-normalized fold change from all IV cohorts as ranks. This approach revealed that communities displayed congruent expression patterns between different cohorts in equivalent comparisons. For example, genes from CM1, CM3, and CM6 were highly induced on day 1 postvaccination, whereas those from CM7 and CM8 were induced on day 7 postvaccination (Fig. 4C).

• Download figure
• Open in new tab
• Download powerpoint

Fig. 4.

Influenza vaccination consensus network. (A) The consensus network was constructed by intersecting CEMiTool results from all IV cohorts and prioritizing frequently detected edges. Inference of communities was performed using a spin glass clustering algorithm. Graph colors are based on community assignment. Each community is represented by a rectangle containing its name and the number of genes and lncRNAs (in parentheses). (B) Overrepresentation analysis of selected communities using blood transcriptional modules (BTMs). False discovery rates (FDR) (x-axis) are plotted for each BTM. (C) Gene set enrichment analysis (GSEA) performed with network communities (rows) and mean fold changes of all vaccinees from each cohort as ranks (columns). Heat map represents normalized enrichment scores (NESs) of communities whose FDR < 0.05. (D) Some PRKCQ-AS1 connections within CM5. The colors of nodes represent log₂ fold-change summaries on day 1 relative to baseline.

Network communities contained dozens of lncRNAs (SI Appendix, Table S2). Meticulous examination of genes connecting to these lncRNAs may unveil potential transregulatory relationships (see SI Appendix). For instance, lncRNA PRKCQ-AS1 is a member of the T cell–related community CM5 in the IV network (Fig. 4A). GSEA revealed that CM5 is consistently down-regulated on day 1 after IV in all IV cohorts (Fig. 4C). A subnetwork of genes connected to PRKCQ-AS1 is also broadly down-regulated on day 1 (Fig. 4D). Several members of this subnetwork are well-known T cell–related genes, such as CCR7, CD27, ITK, LEF1, THEMIS, and SKAP1. These findings suggest that PRKCQ-AS1 may be involved with T cell functions.

The same modular analysis was performed with YF-17D cohorts, which resulted in 17 network communities (SI Appendix, Fig. S6). Several lncRNAs were found within communities related to “chemokines in myeloid cells,” “T cell differentiation,” and “platelet activation” (SI Appendix, Fig. S6). GSEA revealed that, while activation of IFN response was restricted to days 1 and 3 post-IV (CM6 in Fig. 4C), a strong activation of IFN response was still observed on day 7 in the YF network (CM2 in SI Appendix, Fig. S6C).

Differential Expression of lncRNAs in RNA-Seq Data Generated from a New Pediatric Vaccine Cohort.

As RNA-seq analyzes the whole transcriptome, it has the ability to discover lncRNAs associated with vaccination. We assessed whether lncRNAs were differentially expressed in an RNA-seq dataset of 26 children aged 2–4 y old immunized with intranasal LAIV. Vaccination induced differential expression of 55 lncRNAs in whole blood at day 2 post-LAIV compared with baseline, with a false discovery rate (FDR) threshold of 0.05 (SI Appendix, Fig. S7). Several lncRNAs found in our meta-analyses were also found among differentially expressed genes in this cohort, including FAM225A, LRRC75A-AS1, MAPKAPK5-AS1, DANCR, and DICER-AS1.

VaccineDB: A Systems Vaccinology Database.

We developed a user-friendly online database named VaccineDB (https://vaccinedb.sysbio.tools/) to catalog and organize the information generated in this work. By selecting the vaccine (IV or YF-17D), the type of analysis (differential expression, gene-antibody correlation, gene-lncRNA correlation), and gene symbol identifiers, the database displays one forest plot for each time point after vaccination (SI Appendix, Fig. S8). The differential expression and gene-antibody correlation analyses show, respectively, the expression changes and the antibody correlation values of user-input protein-coding genes or lncRNAs after either IV or YF-17D (SI Appendix, Fig. S8). The database also displays the correlation between the expression of a user-provided protein-coding gene and lncRNAs within the same genomic region. This approach may help identify lncRNAs with coordinated expression and potential regulatory cis-acting roles in response to vaccination. Finally, to verify if the lncRNA is associated with known SNPs, the database displays the genomic locus of the lncRNA and all SNPs from the dbSNP database located there.