Virus-specific and shared gene expression signatures in immune cells after vaccination in response to influenza and vaccinia stimulation

Background In the vaccine era, individuals receive multiple vaccines in their lifetime. Host gene expression in response to antigenic stimulation is usually virus-specific; however, identifying shared pathways of host response across a wide spectrum of vaccine pathogens can shed light on the molecular mechanisms/components which can be targeted for the development of broad/universal therapeutics and vaccines. Method We isolated PBMCs, monocytes, B cells, and CD8+ T cells from the peripheral blood of healthy donors, who received both seasonal influenza vaccine (within <1 year) and smallpox vaccine (within 1 - 4 years). Each of the purified cell populations was stimulated with either influenza virus or vaccinia virus. Differentially expressed genes (DEGs) relative to unstimulated controls were identified for each in vitro viral infection, as well as for both viral infections (shared DEGs). Pathway enrichment analysis was performed to associate identified DEGs with KEGG/biological pathways. Results We identified 2,906, 3,888, 681, and 446 DEGs in PBMCs, monocytes, B cells, and CD8+ T cells, respectively, in response to influenza stimulation. Meanwhile, 97, 120, 20, and 10 DEGs were identified as gene signatures in PBMCs, monocytes, B cells, and CD8+ T cells, respectively, upon vaccinia stimulation. The majority of DEGs identified in PBMCs were also found in monocytes after either viral stimulation. Of the virus-specific DEGs, 55, 63, and 9 DEGs occurred in common in PBMCs, monocytes, and B cells, respectively, while no DEGs were shared in infected CD8+ T cells after influenza and vaccinia. Gene set enrichment analysis demonstrated that these shared DEGs were over-represented in innate signaling pathways, including cytokine-cytokine receptor interaction, viral protein interaction with cytokine and cytokine receptor, Toll-like receptor signaling, RIG-I-like receptor signaling pathways, cytosolic DNA-sensing pathways, and natural killer cell mediated cytotoxicity. Conclusion Our results provide insights into virus-host interactions in different immune cells, as well as host defense mechanisms against viral stimulation. Our data also highlights the role of monocytes as a major cell population driving gene expression in ex vivo PBMCs in response to viral stimulation. The immune response signaling pathways identified in this study may provide specific targets for the development of novel virus-specific therapeutics and improved vaccines for vaccinia and influenza. Although influenza and vaccinia viruses have been selected in this study as pathogen models, this approach could be applicable to other pathogens.


Introduction
Although each pathogen has its unique characteristics, pathogen-associated molecular patterns (PAMPs) are conserved across a broad spectrum of pathogens (1). It is well-known that the host relies on limited sets of pattern-recognition receptors (PRRs) to recognize PAMPs, subsequently triggering a variety of signaling pathways as part of the rapid innate immune response (2)(3)(4). Therefore, host responses to pathogens are pathogen-specific although certain immune response pathways are shared across multiple pathogens due to the conservation of PAMPs within those pathogens. While understanding pathogen-specific host responses is critical for the development of pathogen-specific therapeutics (5), identifying shared pathways across a wide spectrum of pathogens can shed light on the molecular mechanisms and/or components which can be targeted for the development of broad or even universal therapeutics and vaccines.
The analysis of gene expression signatures in blood leukocytes has been intensively utilized as a robust approach to characterize host responses after infection or vaccination (6). As such, gene transcripts have been used to distinguish a variety of viral infections (7)(8)(9)(10)(11)(12) or to differentiate viral from bacterial infections (13)(14)(15)(16). This approach has been also applied to study immune responses (17)(18)(19) after vaccination. We have recently applied this approach to compare T-cell transcriptional responses to high-dose and adjuvanted influenza vaccines in older adults (20). A common feature in these studies is that gene expression was characterized in peripheral blood mononuclear cells (PBMCs) after a single pathogen stimulation. Therefore, information on shared pathways of host immune responses to multiple vaccine pathogens is limited and subject to confounding due to differences in experimental/ analytical approaches between studies.
In the vaccine era, almost everyone receives numerous vaccines in their lifetime, often including 1-2 vaccines annually. We surmised that it would be possible to identify shared biological pathways associated with gene expression across multiple pathogens. We also hypothesized that different cell types of blood leukocytes bear unique gene expression signatures. To test this hypothesis, we analyzed gene expression in PBMCs, monocytes, B cells, and T cells isolated from peripheral blood of donors who previously received both influenza vaccine (within <1 year) and smallpox vaccine (within 1 -4 years). We then compared and contrasted vaccine-specific gene expression signatures to identify shared pathways associated with transcriptional responses to the two vaccine pathogens across different cell subtypes.

Materials and methods
The following methods are similar or identical to our previously published studies (20-22).

Ethics statement
This study involved human participants and was approved by The Mayo Clinic Institutional Review Board (IRB# 11-000576). Written informed consents was obtained from each donor before blood collection.

Blood samples and cell separation
Blood samples were obtained from healthy blood donors (n = 10) who received both seasonal influenza vaccine (within <1 year) and smallpox vaccine (within 1 -4 years). Age and gender characteristics of blood donors are summarized in Supplemental  Table S1. The peripheral blood mononuclear cells (PBMCs) were isolated from blood, following our previously reported protocol (22). Isolated PBMCs were suspended in a freezing media (20% heat-inactivated FCS, 10% DMSO in RPMI 1640 media) and stored in liquid nitrogen for future use.

mRNA sequencing
Cells (1×10 6 ) of each population (PBMCs, purified monocytes, B cells, CD8 + T cells) were incubated for 6 hours in a 37°C, 5% CO 2 incubator under three conditions: i) cell culture media only (RPMI 1640 with glutamine supplemented with 10% FCS, Pen/Strep, nonessential amino acid, HEPES buffer) as unstimulated control, ii) influenza virus (A/H1N1/California/07/2009-like strain) in cell culture media with a virus/cell incubation ratio of 5, and iii) vaccinia virus (New York City Board of Health strain) in cell culture media with the virus/cell incubation ratio of 5. After incubation, the cells were collected and subjected to total RNA extraction using a Qiagen RNeasy Plus Mini Kit (catalog no. 74134) and Qiagen RNAProtect reagent (catalog no. 76104), following a step-by-step instruction in the manufacturer's protocol. Extracted RNA was qualitatively and quantitatively evaluated by Agilent 2010 Bioanalyzer assay (Agilent, Palo Alto, CA) and a NanoDrop 2000 spectrophotometer (Thermo Fisher Scientific, Waltham, MA).
The extracted RNA (150 ng) was sequenced at Mayo Clinic's Gene Sequencing Core Facility (Rochester, MN). Briefly, the TruSeq ® Stranded mRNA Library Prep v2 (Illumina, San Diego, CA, USA) was used to create cDNA libraries, following the manufacturer's protocol used at Mayo Clinic Core Facility. Illumina's HiSeq 2000 S2 Reagent Kit (51 cycles) was used to perform paired-end read on the Illumina HiSeq 2000 Instrument. Gene sequencing data was aligned using the MAP-RSeq V1 pipeline to the h19 human genome as outlined in Kalari et al. (23). The viral sequences were mapped to the vaccinia virus ACAM2000 genome (AY313847.1) and the influenza A/ California/07/2009 H1N1 genome.

Statistical analysis
Gene expression was normalized for each cell subset (PBMCs, purified monocytes, B cells, CD8 + T cells) using Conditional Quantile Regression (CQN), which adjusts for GC content and gene length bias, via the cqn R package (24). Any post-normalized genes with <16 counts and genes with coefficient of variation in the lower 20 th percentile in all groups were removed from further analysis. After filtering, 12,733, 11,641, 12,201, and 12,501 human genes remained for the PBMCs, monocytes, B cells, and T cells, respectively. The edgeR package (25,26) was used to identify differentially expressed genes (DEG) by fitting quasi-likelihood negative binomial generalized log-linear model, utilizing the offset provided by the CQN normalization and blocking on subject. DEGs were identified based on fold change (FC) of gene expression as compared to corresponding uninfected samples with a threshold set at abs(log2FC) ≥0.585 and p-value <0.05. All p-values associated with DEG results were unadjusted. Gene set enrichment analysis (GSEA) and over-representation analysis using KEGG pathways was performed using the clusterProfiler package in R (27). Pathway analyses were performed under the settings: i) the minimum geneset size of 3, ii) the maximum geneset size of 800, iii) the pvalue cutoff of 0.05, and iv) adjusted p-values were calculated using the Benjamini-Hochberg method (28). GSEA provides a statistical framework using the hypergeometric distribution to test if KEGG pathways are overrepresented by significant DEGs.

Differentially expressed genes in influenza-stimulated cells
After blood collection, we isolated PBMCs, and purified monocytes, B cells, and CD8 + T cells were using commercial magnet-based separation kits ( Figure 1). We incubated each of these cell populations with either influenza virus or vaccinia virus and profiled gene expression in each cell subset ( Figure 1). Differential gene expression analysis identified a total of 2,906, 3,888, 681, and 446 differentially expressed genes (DEGs) in influenza-stimulated PBMCs, monocytes, B cells, and CD8 + T cells, respectively (Figures 2A-D). Interestingly, PBMCs and monocytes shared 1,763 DEGs, while PBMCs had 322 and 284 DEGs in common with B cells and CD8 + T cells, respectively ( Figure 2E). In overlapping DEGs, we observed a similar pattern of gene expression in PBMCs and monocytes ( Figure 2F). The expression patterns of DEGs were also similar between PBMCs and B cells ( Figure 2G), and between PBMCs and CD8 + T cells after influenza stimulation ( Figure 2H), although some genes differentially regulated in B cells and CD8 + T cells.
Among the DEGs in each cell population, the highest fold changes of gene expression were observed in PBMCs and monocytes (Table 1). PBMCs and monocytes also shared numerous DEGs, including seven of top ten up-regulated genes (IFNA1, IFNA2, IFNA7, IFNA10, IFNA13, IFN17, TEKT1) and two of top ten down-regulated genes (CCL24 and SCRT2) ( Table 1). In contrast, the expression level of DEGs in B cells and T cells werẽ 100-fold lower than that in PBMCs and monocytes (Table 1). Both B cells and T cells did not share any DEGs with PBMCs among their top DEGs (Table 1). Meanwhile, B cells and T cells shared some genes encoding for interferon induced proteins, such as IFIT1, RSAD2, among their top DEGs (Table 1).
To further understand the biological pathways associated with the DEGs, we performed gene set/pathway enrichment analysis using the clusterProfiler package and KEGG modules and pathways (27). We found that DEGs in PBMCs, monocytes, B cells, and CD8 + T cells were significantly enriched in 42, 30, 11, and 26 pathways (adjusted p <0.05), respectively ( Figure 3). Influenza A and coronavirus disease-COVID-19 pathways were among the top enriched pathways in all influenza-stimulated cell populations ( Figure 3). The NOD-like receptor signaling pathway was also enriched in each of the four cell populations upon influenza stimulation ( Figure 3). Similar to the pattern of gene expression ( Figure 2F), almost all signaling pathways associated with DEGs in monocytes were also found in PBMCs, including JAK-STAT signaling pathway, cytosolic DNA-sensing pathway, Toll-like receptor signaling pathway, RIG-I-like signaling pathway, cytokine-cytokine receptor pathway, and PI3K-Akt signaling pathway ( Figure 3).

Differentially expressed genes in vaccinia-stimulated cells
We also identified DEGs in vaccinia-stimulated cells ( Figure 4). In general, the number of DEGs in vaccinia-stimulated cells was approximately 30-fold less than that in the influenza-stimulated cells. Specifically, a total of 97, 120, 20, and 10 DEGs were identified in PBMCs, monocytes, B cells, and CD8 + T cells after vaccinia stimulation, respectively ( Figures 4A-D). Of vaccinia-induced DEGs, 36 DEGs were identified in both PBMCs and monocytes ( Figure 4E). We also observed a similar expression pattern of overlapping DEGs in PBMCs and monocytes ( Figure 4F), PBMCs and B cells ( Figure 4G). Only one common DEG was identified in PBMCs and CD8 + T cells, but its expression was in opposite directions ( Figure 4H). Similar to influenza stimulation (Table 1), vaccinia-stimulated PBMCs and monocytes shared multiple interferon-encoding genes, such as IFNA7, IFNA17, among their top DEGs (Table 2).
Gene set enrichment analysis demonstrated that vacciniainduced DEGs in PBMCs, monocytes, and CD8 + T cells were significantly (adjusted p <0.05) enriched in 64, 29, and 9 pathways ( Figure 5). No pathway was significantly associated with vaccinia-induced DEGs in B cells. Similar to influenza-induced DEGs ( Figures 3A, B), vaccinia-induced DEGs in PBMCs and monocytes were enriched in multiple innate pathways, including cytosolic DNA-sensing pathway, RIG-I receptor signaling pathway, Toll-like receptor signaling pathway, JAK-STAT signaling pathway, and NOD-like receptor signaling pathway ( Figures 5A, B). Study design. Healthy subjects (n = 10) who received both seasonal influenza vaccine (within the prior year) and smallpox vaccine (with the prior 1-4 years) were enrolled for this study. Blood was sampled from each subject and PBMCs were isolated from blood. Monocytes, B cells, and CD8 + T cells were separated from PBMCs using MACS MicroBead isolation kits. Each cell population was then incubated with either influenza virus (A/H1N1/ California/07/2009-like strain) or vaccinia virus (NYCBOH strain). Cells in cell culture media served as unstimulated controls. After infection, RNA was extracted from incubated cells and gene expression in each cell population was profiled.

Shared DEGs in each cell population after influenza and vaccinia stimulation
We next explored the commonalities of DEGs in each cell population after influenza and vaccinia stimulation. Of DEGs induced either by influenza or vaccinia virus, we identified a total of 55, 63, and 9 shared DEGs in PBMCs, monocytes, and B cells, respectively ( Figure 6). Meanwhile, no shared DEGs were identified in CD8 + T cells. Although we observed a similar expression pattern in PBMCs and monocytes after viral stimulation; a significant number of shared DEGs were expressed in different directions upon influenza and vaccinia stimulation (Figures 6A, B). In contrast, shared DEGs in B cells were expressed similarly after influenza and vaccinia stimulation ( Figure 6C).
We further identified over-represented pathways associated with the shared DEGs after influenza and vaccinia stimulation to understand their biological functions. We found that shared DEGs in PBMCs and monocytes were significantly over-represented in 35 and 26 biological pathways, respectively (Figure 7). These pathways included: cytokine-cytokine receptor interaction, Toll-like receptor signaling pathway, NOD-like receptor signaling pathway, RIG-I-like receptor signaling pathway, viral protein interaction with cytokine and A B D C FIGURE 3 KEGG pathways associated to DEGs in influenza-stimulated PBMCs (A), monocytes (B), B cells (C), and CD8 + T cells (D). All KEGG pathways shown are statistically significant (adjusted p <0.05). NES, normalized enrichment score. Red color (NES <0) indicates suppressed pathways while bule color (NES >0) indicates activated pathways. The size of each red or blue dot represents the total number of genes (size of geneset) in the pathways. Additional information of these pathways is provided in Supplemental Tables S10-S13.

Discussion
Understanding host transcriptomic responses to pathogens after vaccination is critical for the development of vaccines (5). In this study, we identified virus-specific DEGs in PBMCs, monocytes, B cells, and T cells isolated from healthy donors who had been immunized with both seasonal influenza vaccine and smallpox vaccine. Among the virusspecific DEGs, we were able to identify shared DEGs and their associated biological pathways after influenza and vaccinia virus stimulation. Our data showed a significant number of innate signaling pathways associated with shared DEGs were activated in response to influenza and vaccinia virus stimulation. Our data also highlighted the role of monocytes as a cell population driving transcriptomic responses in PBMCs to both viruses. These results suggest that this specific cell type may be useful as sentinels to monitor the physiologic response during an infection or after vaccination.
We chose influenza and smallpox vaccines as vaccine models for multiple reasons. First, influenza vaccine is administered annually in the U.S.; therefore, it is a logical choice to study. Meanwhile, after its eradication in 1980, routine smallpox vaccination ceased (29). However, vaccinia virus is used as a vector for the development of multiple vaccines. Public concern still exists concerning the safety of these vaccines and there have been an increased number of zoonotic orthopoxvirus outbreaks. A recent global re-emergence of monkeypox (30), another member of Poxviridae, further heightens this public concern. As a result, vaccinia virus serves as a re-emerging pathogen model which may require a vaccination campaign to re-establish herd   Supplemental Tables S6-S9. immunity (31,32). Second, by choosing two distantly related viruses, such as influenza (a small RNA virus that replicates in the nucleus) and vaccinia (a large DNA virus that replicates in the cytoplasm), we aimed to demonstrate that shared signaling pathways of host responses exist across highly diverse pathogens.
In other words, we demonstrated that some shared signaling pathways were activated in response to the stimulation of a wide spectrum of pathogens. Therefore, these signaling pathways should be principal pathways in response to viral infections and identifying these pathways will aid the development of general or even universal therapeutics and vaccines.
Previously, we profiled transcriptomic responses in blood leukocytes after influenza vaccination (21, 33-37) and smallpox vaccination (38-42). We and others repeatedly observed the activation of interferon-encoding genes in PBMCs after either influenza (35-37, 43) or vaccinia stimulation (38, 40). Consistent to these previous findings, in this study we found that interferonencoding genes were up-regulated in both PBMCs and monocytes stimulated with either influenza virus (Table 1) or vaccinia virus  (Table 2). Transcriptional analysis has also found that interferonencoding genes were among the top activated genes in human PBMCs infected with different viruses (7,12). Together, these findings A B C FIGURE 5 KEGG pathways associated to DEGs in vaccinia-stimulated PBMCs (A), monocytes (B), and CD8 + T cells (C). All KEGG pathways shown are statistically significant (adjusted p <0.05). NES, normalized enrichment score. Red color indicates suppressed pathways (NES <0) while bule color indicates activated pathways (NES >0). The size of each red or blue dot represents the total number of genes (size of geneset) in the pathways. Additional information of these pathways is provided in Supplemental Tables S14-S16. consolidate a defensive role of interferons in viral infections. Serving as the first line of innate defense, interferon production is amongst the first innate responses against viral infection (44, 45); hence, this could explain for the up-regulation of interferon-encoding genes observed in this and other studies. In turn, the production of interferons is activated by the JAK-STAT signaling pathway (44-46), which was consistently amongst the top enriched signaling pathways in this study (Figures 3, 5,  7). Altogether, these findings suggest the JAK-STAT signaling pathway as a main target for the development of viral therapeutics or vaccines (46). After influenza stimulation, the identified DEGs were enriched in influenza A and coronavirus disease-COVID-19 pathways in each cell population ( Figure 3). These results suggest a similarity between these two respiratory viruses in the way they are recognized and targeted by the host immune system, as observed by others across respiratory viruses (9, 47). More importantly, even with the influenza and vaccinia viruses used in this study, we could identify 35 and 26 biological pathways associated with shared DEGs in PBMCs and in monocytes, respectively (Figure 7). These shared pathways were mainly involved in the activation of innate immune responses, such as NOD-like, Toll-like, RIG-I-like receptor signaling pathways, cytokine-cytokine receptor interaction, cytosolic DNA-sensing pathway, and natural killer cell mediated cytotoxicity (Figure 7). Altogether, these observations suggest that the host relies on a similar set of PRRs to sense the presence of pathogens and trigger similar innate immune pathways in responses to pathogens (1)(2)(3)(4)48). Therefore, these innate immune pathways could potentially be targeted for the development of broader viral therapeutics or be used to inform the use of appropriate adjuvants for vaccines.

A B C
Consistently, we observed the similarity in DEGs and their expression patterns in PBMCs and in monocytes upon either influenza or vaccinia (Tables 1, 2; Figures 2A, B, F, 4A, B, F, 6A, B). As a result, the biological pathways associated with the DEGs in PBMCs and in monocytes were also similar ( Figures 3A, B, 5A, B, 7A, B). These results reflect the importance of monocytes as a cell population driving transcriptomic responses. This may be due to the remarkable degree of plasticity observed in myocytes, in which they can rapidly sense and respond to diverse invaders (49, 50). Our results suggest that monocytes could be a major target cell type for viral therapeutics and/or serve as a cellular biomarker for infection/ vaccination. In fact, previous studies have found that monocytes are the main target for a variety of viral infections, including influenza, dengue, measles, Zika (51)(52)(53)(54). Our results further demonstrate the heterogenous nature of PBMC samples and the limitations in using whole blood or PBMCs to evaluate transcriptomic changes after infection or vaccination. While useful as an initial step, examination of isolated cell subsets or the use of single cell technologies are more powerful and increasingly taking the place of bulk RNA-Seq.
In contrast, B cells and CD8 + T cells provided limited information on transcriptomic responses and the biological pathways associated with DEGs in B cells and CD8 + T cells did not show a specific trend (Figures 2-5). This could be due to the fact that receptors on B cells and T cells are binding to specific viral epitopes, rather than DAMPs. B cells and T cells are involved in adaptive immune responses, which chronologically occurs after the activation of innate responses triggered by PRRs.
Our study has both strengths and limitations. The strengths of this study mainly lie in our study design. As such, first we studied the transcriptomic responses induced by two viral vaccine pathogens, A B FIGURE 7 KEGG over-representation pathways associated to overlapped DEGs in PBMCs (A) and monocytes (B) after influenza and vaccinia stimulation. All KEGG over-presentation pathways shown are statistically significant (adjusted p <0.05). The size of each dot represents the count of DEGs in associated pathways. The adjusted p value (p.adjust as shown in the figures) is color coded. Additional information of these pathways is provided in Supplemental Tables S17-S19.
which minimizes variation and confounding due to the differences in experimental/analytical approaches between studies. To the best of our knowledge, this study is one of the first studies identifying shared transcriptomic responses induced by these two vaccine pathogens in the same host. Second, rather than study only PBMCs we expanded the transcriptomic characterizations in purified monocytes, B cells, and T cells. By studying the transcriptomic responses in these purified cell subset populations and comparing the responses in these cell populations with the responses in PBMCs, we found that the transcriptomic responses in PBMCs were driven by monocytes. Meanwhile, our observations were limited to 10 subjects, mainly due to the difficulty in recruiting study participants who met the requirements for this study. A future study with the same design, but with a larger study cohort may consolidate the results observed in this study. In this study, we used only two virus models. A future study including more pathogens of different types, such as bacteria, virus, fungi would facilitate the development of universal therapeutics or vaccines for a broader spectrum of pathogens. In addition, we used fold change of gene expression and unadjusted p value to identify DEGs. While the fold change shows a relative expression of each individual gene in virus-infected cells as compared to control cells, unadjusted p value may be impacted by multiple testing, where the expression of a certain gene is influenced by others (55).
In conclusion, we identified both virus-specific and shared DEGs in response to influenza and vaccinia stimulation and their associated biological pathways in this study. Our results highlight monocytes as a main cell type driving transcriptomic responses in the PBMCs. The shared biological signaling pathways observed in this study may provide molecular insights into virus-host interactions, serving as targets for the development of novel virusspecific therapeutics and new generic vaccines for vaccinia and influenza. Although influenza and vaccinia viruses have been selected in this study as pathogen models, we believe that this approach is applicable to other pathogens including SARS-CoV-2.

Data availability statement
The data presented in the study are deposited in Synapse, under Project SynID syn52165181.

Ethics statement
This study was approved by The Mayo Clinic Institutional Review Board. Each study participant provided their written informed consents before blood collection.

Author contributions
HQ analyzed data, wrote, and revised the manuscript. KG and DG analyzed mRNA sequencing data and reviewed the manuscript. RK, IH, and IO carried out the experiments, and critically reviewed and edited the manuscript. GP and RK acquired research funding, supervised the project, and edited the manuscript. All authors contributed to the article and approved the submitted version.

Funding
This work was supported by the NIH through the NIAID Population Genetics Analysis Program Contract No. HHSN266200400065C and Contract No. HHSN272201000025C, and by the National Institutes of Health Grants U01 AI89859 and R01 AI132348. The funding agencies had no role in study design, data collection and analysis, decision to publish, or preparation of this manuscript. The remaining authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher's note
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.