|Home | About | Journals | Submit | Contact Us | Français|
To identify sources of inter-subject variation in vaccine responses, we performed high-frequency sampling of human peripheral blood cells post-vaccination, followed by a novel systems biology analysis. Functional principal component analysis was used to examine time varying B cell vaccine responses. In subjects vaccinated within the previous three years, 90% of transcriptome variation was explained by a single subject-specific mathematical pattern. Within individual vaccine response patterns, a common subset of 742 genes was strongly correlated with migrating plasma cells. Of these, 366 genes were associated with human plasmablasts differentiating in vitro. Additionally, subject-specific temporal transcriptome patterns in peripheral blood mononuclear cells identified migration of myeloid/dendritic cell lineage cells one day after vaccination. Upstream analyses of transcriptome changes suggested both shared and subject-specific transcription groups underlying larger patterns. With robust statistical methods, time-varying response characteristics of individual subjects were effectively captured along with a shared plasma cell gene signature.
Pandemic and seasonal influenza causes significant morbidity and mortality worldwide. A key strategy in preventing influenza infection is vaccination, which induces B cells to differentiate through several stages into plasmablasts and long-lived plasma cells. These cells produce protective antibodies (immunoglobulin) against hemagglutinin (HA) capsid proteins1. The World Health Organization (WHO) defines effective influenza vaccine responses as 28-day post-vaccine antibody hemagglutinin inhibition assay (HAI) titer ≥ 1:402. However, influenza vaccines produce highly variable B cell responses among individuals3, making it difficult to predict who will achieve protective antibody titers and when.
Advances in data-driven personalized medicine, as applied to vaccine immunobiology, require new statistical methods that are able to capture immune response features common across subjects (e.g. consensus gene expression patterns) and patient-specific patterns which differ between subjects. This is especially true in vaccine studies. For example, recent studies on transcriptome patterns in vaccine responses have used sparse sampling (i.e. microarray analysis of PBMC on days 0, 3, and 7), and analytic methods that do not involve time series analysis4,5,6. In this report, we describe a robust systems biology approach using functional principal component analysis to identify common and subject-specific time-varying gene signatures within the trivalent influenza vaccine (TIV) immune response.
In a memory response to influenza vaccination, resting memory B cells differentiate through proliferating plasmablast states to terminally differentiated non-proliferating plasma cells7. Migrating through peripheral blood to bone marrow and spleen, virus-specific antibody-secreting cells (ASC) appear 5–9 days post-vaccination8. Vaccine or influenza exposure induces memory B cells producing cross-reactive antibodies to similar influenza strains, termed heterosubtypic immunity9. The presence of heterosubtypic immunity alters the speed and strength of TIV responses. Previous systems biology vaccine response analyses have sampled subjects on a single day empirically described as the peak B cell response10,4. Such approaches lack the structure to capture peak responses or the statistical power to separate time varying gene signatures from inter-subject variation. Our approach is fundamentally different. Our goal was instead to identify differences in the character of the transcriptome responses of individuals using high-frequency sampling and an analytical method designed to characterize temporal patterns in complex systems. In contrast to other studies10,4, we were not trying to predict who responded to the vaccine. Rather, our goal was a deeper understanding of the time varying, vaccine specific, transcriptome response, and to determine if the B and plasma cell gene expression patterns varied as a function of pre-existing vaccine antigen memory. Systems analysis of high-frequency time course data requires more samples per subject, but better defines peak and time varying B cell immune responses.
To separate concensus from subject-specific time-varying B cell transcriptome responses to influenza vaccination, we performed high-frequency (daily) sampling of serum, peripheral blood mononuclear cells (PBMC), B cells and plasma cells from 14 human subjects over 11 days post-TIV administration. Our approach combined time-series transcriptome analysis, flow cytometry, ELISPOT, ELISA, and HAI assays with computational analysis to identify time-varying RNA biomarkers associated with robust B cell responses. To identify subject-specific gene expression signatures, we used functional principal component analysis (FPCA), a statistical dimension reduction method for high frequency time series data similar to principal component analysis for single time-point data11,12. In essence, the differential expression of each gene is determined by its own signal-to-noise ratio. The major temporal variation in the RNAseq time-series data can be characterized by a few representative patterns called eigenfunctions. Each gene expression curve is approximated by a linear combination of these eigenfunctions, enabling a flexible and parsimonious representation of time-series data. This method accommodated individuality in vaccine response by identifying time-varying transcriptome patterns unique to each subject. For a detailed explication of FPCA, and underlying statistical theory, please see the specific section in Supplementary Methods.
This approach identified a peripheral blood B cell transcriptome signature in previously vaccinated subjects that was strongly correlated with migrating plasma cells observed with cytometry. Upstream analysis of significant gene sets revealed common immunological dynamics between subjects as well as unique vaccine response features in individuals. At the subject level, these responses were associated with robust protective antibody responses and a prior history of influenza vaccination. Eigenfunctions calculated for PBMC gene expression also demonstrated an early wave of myeloid/dendritic cell lineage migration peaking 24 hours post-vaccination. Additional eigenfunctions identified transcriptome changes unique to each individual, commonalities between subject-specific signatures, and their relative contribution to the overall vaccine response.
We vaccinated three cohorts of 3–6 subjects (n = 14 subjects) with the 2010–2011 TIV. Subjects reported mixed vaccine histories, with two subjects reporting no influenza vaccination in three consecutive years prior to the study, and three subjects reporting receiving all available influenza vaccines. Two of study vaccine strains, A/California/7/2009 and B/Brisbane/60/2008 were included in recent vaccines (Supplementary Table S1).
For analysis, peripheral blood was drawn during the week prior to vaccination (pre-V), immediately before vaccination (day 0), daily for days 1–10 and on day 21 post-vaccination. HAI and isotype-specific ELISA assays were performed on serum samples. (Fig. 1, Supplementary Fig. S1). No two subjects had identical changes in HAI titers. One subject had no pre-vaccine immunity (HAI < 1:40), and five subjects had pre-vaccine immunity (HAI ≥ 1:40) to all three viral strains (Fig. 1, Supplementary Table S2). At day 21, all 14 subjects had immunity to A/California, 12 had immunity to A/Perth and 13 had immunity to B/Brisbane.
B cell samples (days 0–10) were analyzed by RNAseq, flow cytometric and ELISPOT (880 samples) (Supplementary Fig. S1). Due to cost, RNA Seq analysis was performed on PBMC and B cell enriched samples from one contemporaneous cohort of five subjects, days 0 to 10 (110 samples). S01 RNA was used for RNA Seq pilot studies. We selected a single cohort of subject samples for RNA Seq analysis. Each gene was tested for differential change from day 0 to day 10 using FPCA-based significance testing13 and a permutation test coupled with multiple-testing correction to control the false discovery rate (FDR < 0.05)14. A deviation in a gene expression curve from a flat line indicates a significant temporal pattern in transcriptomic variation.
Intra-subject variation in differentially expressed genes was moderate (p-values in Supplementary Table S3). S02 had the largest number of significantly time-varying genes (5580 genes) while S06 had few (2 genes). For better comparison with other subjects, 1600 genes with the highest signal to noise ratio were added to the analysis of S06. After excluding genes with less than a 1.2-fold-change from baseline, the numbers of temporally differentially expressed genes still varied greatly: S02 (5256 genes), S03 (1309 genes), S04 (2147 genes), S05 (1603 genes) and S06 (1053 genes). Notably, 742 genes were common to three of five subjects.
We next used FPCA to identify the major variations in temporal RNA expression patterns of the differentially expressed genes for each subject, which were then represented by eigenfunctions (Fig. 2a, complete list of significantly time-varying genes in Supplementary Table S3). For three subjects (S02, S03, S04), 90% of transcriptome variation could be accounted for by a single eigenfunction. All three of these subjects had a previous history of influenza vaccination, pre-vaccine immunity to at least one viral vaccine strain by WHO criteria, and strong antibody responses by HAI titer (Fig. 2b, p-values in Supplementary Table S2) and vaccine-specific ELISA (Fig. 2c, p-values in Supplementary Table S2).
Cellular responses were seen by FACS in the appearance of vaccine-specific ASC (Fig. 2d), and migrating CD27hiCD38hiCD20lo CD138+ plasma cells and CD27hiCD38hiCD20loCD138− plasmablasts (Fig. 2e, p-values in Supplementary Table S4). In semi-quantitative ELISPOT data, (Supplementary Fig. S2), we found secretion-rate profiles typical of activated ASC15, with high-rate vaccine-specific IgG-secreting cells appearing at peak responses.
Subject-specific dynamics were seen in vaccine response characteristics at multiple levels. S02 had a sharp peak in cellular and RNA changes at day 5, whereas S03 and S04 had more gradual changes. These subjects met serum and cellular criteria for robust memory B cell responses to TIV. However, with some evidence of de novo reactivity to A/California, and without specific clone data, we prefer the term “predominantly recall” responses.
In two subjects (S05, S06), 90% of the variation in gene expression could be explained by two eigenfunctions, indicating two distinct temporal patterns of gene expression. S06 had a linear trend in some genes and S05 had multimodal eigenfunctions. These subjects reported no influenza vaccination in the previous three years and developed de novo immunity to at least two vaccine-specific viruses by WHO HAI criteria. Acknowledging the difficulties in defining an influenza response in adults as truly naive, we prefer the term “predominantly de novo” responses.
Taken together, results from the high frequency sampling data show coordinated and dynamic vaccine-response changes at serum, cellular, and transcriptomic levels, as well as inter-subject variations in immunological response characteristics. Features of time-varying transcriptome profiles independently identified subjects with significant pre-vaccine B cell memory. Also, high-frequency sampling provided the resolution to discern unique expression dynamics across assays for each subject.
To connect the transcriptomic and cellular level B cell influenza vaccine responses, we hypothesized that predominantly recall TIV responses would be characterized by significant variations in established markers of late B cell differentiation (Fig. 2f). Detailed analysis showed known genes with time-varying patterns that paralleled B cell phenotypic changes in S02, S03 and S04 including the genes for CD70 receptor CD2716,17,18, ecto-enzyme CD3816, plasma cell marker syndecan-1 (SDC1 or CD138), and genes encoding regulatory molecules PRDM1 (Blimp-1)19,7, IRF416,20,21 and XBP122. Also increased was CD59, which has been associated with memory B cell differentiation23. Increased expression of CXCR3, a gene for a migration marker24 suggested trafficking from the germinal center through peripheral blood presumably to bone marrow and spleen. Decreases in MS4A1 (CD20) and CD19 at peak cellular changes are also consistent with late-stage differentiating B cells16,25,26,23. PAX57 and BACH227 are known repressors of B cell differentiation and decreased gene expression has been linked to expression of PRDM128. In subjects with predominantly de novo responses (S05, S06) time-varying expression of these differentiation markers was modest or absent.
Decreases were seen at peak response in CD5, the gene expressing a marker of naïve B cells, and SERPINB9 (PI-9), an enzyme that protects cytolytic B cells from endogenously produced granzyme B29. It has been proposed that a granzyme B expressing cytolytic B cell lineage diverges from plasmablasts during differentiation30. Reduced expression of these markers may reflect ASC fractional expansion post-vaccination.
To confirm and validate changes in expression of key B cell differentiation genes, we probed RNA samples of the eight other subjects by qRT-PCR for 11 genes: SDC1, IGHG, XBP1, CD27, CD38, CD59, PRDM1, IRF4, IRF8, BACH2, and PAX5. These subjects were heterogeneous in vaccine response by HAI titer (Fig. 3a), vaccine-specific ELISA (Fig. 3b), and flow cytometry (Fig. 3c). Changes in migrating B cells and plasma cells over time exhibited a similar pattern in these independent subjects (Fig. 3d and Supplementary Table S5).
In comparing significant time-varying gene sets between individuals, we found a biomarker signature of 742 genes common to predominantly recall subjects (Fig. 4, p-values in Supplementary Table S6). To identify statistically significant connections between transcriptomic and cellular levels, we compared phenotypic and gene expression data using Pearson's correlation test with the Benjamini-Hochberg multiple test correction to control FDR at 0.05. Out of 742 genes shared among the previously vaccinated individuals, 740 (S02), 404 (S03), and 238 (S04) genes were significantly correlated with CD27hiCD38hiCD138− plasmablasts, indicating a very strong association. This gene set was even more highly correlated with CD27hi CD38hiCD138hi plasma cells. Out of 742 shared genes, 742 (S02), 553 (S03), and 725 (S04), were correlated with plasma cells (FDR < 0.05) (Fig. 5a, correlations in Supplementary Table S6). These associations were strong in all three subjects, even though the timing of peak expression, and composition of genes within the broader significant gene sets, varied between subjects. We refer to this set of 742 genes as the plasma cell gene signature (PCgs).
We next asked if the PCgs could be validated by examining transcriptome data from in vitro human plasma cell differentiation. Peripheral blood CD27+ IgG+ memory B cells from six normal healthy volunteers were CFSE labeled to track proliferation and stimulated in vitro for 60 hrs with CpG2006 ODN, IL-2, IL-10, IL-15, and BAFF23. For transcriptome analysis, cells were FACS-sorted into three stages of late-stage memory B cell activation and differentiation; undivided CD27lo memory B cell, proliferating CD27lo cells and proliferating CD27hi plasmablasts (Fig. 5b) (n = 6 subjects). We hypothesized that a subset of PCgs genes would be differentially expressed between the CD27lo (proliferating but non-IgG secreting) B cell and the more differentiated CD27hi plasmablasts (Fig. 5c). Of 2033 genes differentially expressed between CD27lo and CD27hi B cells, 366 genes of the PCgs were highly expressed in the CD27hi plasmablasts (Fig. 5d, full gene list in Supplementary Table S6). These independent findings strengthen the association between the PCgs and differentiating plasma cell subtypes (Fig. 5e).
Given the strength of the in vivo and in vitro associations between the plasma cell population kinetics and the PCgs at the gene level, it seemed likely that there were common upstream regulatory transcription factors. We performed upstream analysis on gene expression data from S02 using Ingenuity® Pathways Analysis (IPA). Fisher's exact test was used to assess overlap between the data set and curated gene sets associated in the literature with upstream modulators (www.ingenuity.com). We constructed a time series movie of a network of identified upstream modulators (p-value < 0.01) and their targets, (see Supplementary Movie S1). Gene expression changes consistent with release of suppression of B cell regulator PRDM1 (Blimp1) by BACH227 occur by day 4. In days 4 through 10, upregulation of XBP1 was seen, as well as downstream genes associated with ER function, the unfolded protein response, protein production and transport. We also performed functional enrichment analyses of Gene Ontology (GO) terms using DAVID (http://david.abcc.ncifcrf.gov/) and analyses of pathways with KEGG, BIOCARTA and REACTOME (Supplementary Table S7). Ten of the top 30 categories of functionally related genes in the PCgs involved ER function and protein production. These findings are consistent with involvement of the PCgs in program-level upregulation of antibody production machinery and the unfolded protein response seen during plasma cell development28.
In addition to the common PCgs, each subject had unique time-varying gene sets. We hypothesized that these may reflect differences between subjects in regulatory transcription factors. Using the complete set of significant genes for each predominantly recall subject, upstream analysis identified a likely set of upstream activators (Z-score > 2.0) that spanned the peak B cell response in all three subjects (Fig. 6). IL2, IL4, IL5, IL6 and CD40 ligand were identified as probable activators. While the signaling pathways of these activators are not mutually exclusive, this is consistent with a plasma cell history of germinal center reaction.
In addition to common upstream activators, each subject also had a unique set of significantly varying transcription factors. Of particular interest, a wave of probable activation of Type I IFN and IFN-related gene sets appeared in S04 two days after the shared germinal center pattern. This could be evidence of a secondary migration of B cells with previous plasmacytoid dendritic cell (pDC) interaction within the broader peak response. Producers of large amounts of Type I Interferons, pDC are part of the innate response to virus and can induce antibody production by CD40-stimulated B cells in response to influenza31,32. As interferon regulatory factor (IRF) proteins can also regulate B cell differentiation in response to toll-like receptor stimulation33,34, another interpretation of this pattern is a viral DNA-initiated innate response to the vaccine. The eigenfunctions derived by the analysis, along with upstream analysis, identified both common and subject-specific dynamic gene sets that may represent the effects of interacting innate and adaptive cellular immune systems on differentiating plasma cell.
We next hypothesized that the PCgs could be identified in PBMC collected in parallel with the B cell samples. We found something quite unexpected. The primary eigenfunction in the RNA Seq gene expression data, which explained over 70% of the variation for four out of five subjects (S02, 81.44%; S03, 86.68%; S04, 72.8%; S05, 57.15%; S06, 85.32%), peaked on day 1 post-vaccination (Fig. 7a). In the predominantly recall subjects, the genes with largest positive loadings on the primary eigenfunction were enriched for myeloid/dendritic cell (DC) lineage-associated genes (Fig. 7b, Supplementary Table S8). Interestingly, the first eigenfunctions for predominantly de novo response subject S06 also contained myeloid/DC genes with more modest overall gene expression changes. These results suggest that even within the far more complex and dynamic mix of cellular populations in post-vaccination PBMC, FPCA analysis could identify a wave of myeloid lineage cells released into peripheral blood one day post-influenza vaccination in four of five subjects.
Using a method based on FPCA for time-series transcriptome analysis, we identified a 742 gene signature, upstream transcription factors, and subject-specific gene sets associated with a robust B cell response to influenza vaccine. The method yields quantitative functions (eigenfunctions) that describe the temporal pattern of major transcriptome variation. In previously vaccinated subjects, the 742 gene signature was highly correlated with the migration of differentiating plasma cells through peripheral blood. There was an association between the PCgs and key upstream transcription factors that control antibody production including PRDM1 and XBP1. A second gene signature was identified in unfractionated PBMC with a transcriptome profile suggesting migration of cells of the myeloid/DC lineage one day post-vaccination. In addition to these common gene signatures, eigenfunctions identified specific points of variation between individuals including a late spike in IFN-related gene activity. Thus, this approach can identify processes common across individuals in a vaccine response, as well as unique features that may underlie variation between subjects.
A systems biology approach to studies of B cell response requires not only connections between nested levels of regulation from genes to the whole organism, but also perturbations in these systems over time. This study was designed to capture frequent and complete B cell gene expression dynamics after influenza vaccination in a small set of subjects. In this study, high-frequency data from each subject serves as an intra-subject control, adding statistical power and accuracy for detection of significant signals from noisy transcriptome data. A better understanding of vaccine-induced subject-specific gene dynamics may help design studies with more subjects and fewer timepoints that still capture transcriptome patterns.
The use of eigenfunctions to identify common patterns of gene expression that explain transcriptome variations is a powerful analytic method. Inter-subject vaccine response heterogeneity is considerable3,9,4. To generalize conclusions, other systems biological studies have averaged gene expression data across subjects10 or restricted their analyses to significant genes shared in at least 20% of subjects4. Our goal was not to predict vaccine efficacy in large groups of people and indeed, all of our subjects achieved protective HAI titers. Our goal was to identify differences between the vaccine responses of individual subjects over time, and to identify dynamic gene expression mechanisms that underlie those differences. Our method approached subject variation by identifying subject-specific time-varying transcriptome signatures and within those unique patterns, identifying a B cell response gene set shared among individuals. This approach may yield more robust insights into immune response dynamics than sampling at wide time intervals (e.g. only days 0, 4 and 7 post-vaccination), which misses critical events in the individual response. While other studies have found plasma cell associated genes upregulated on Day 76,4, we found subject-specific variation in peak expression of PCgs that would have been missed with typical interval sampling.
Using high-frequency sampling, we found that subjects that had been previously vaccinated each had individual response characteristics, but 90% of the gene expression variation over time could be accounted for by a single subject-specific eigenfunction containing a common plasma cell transcriptome signature. This gene signature was identical in all three of the previously vaccinated subjects. As our sample size is small and all of our subjects for RNA Seq achieved protective HAI titers for the vaccine antigens, it is difficult to compare vaccine efficacy with incidence of the PCgs. One interpretation suggested by our findings is that the predominantly recall influenza vaccine responses exhibited a synchrony of B cell differentiation across vaccine-specific memory B cell populations, while the predominantly de novo responses did not. In combination with stochastic nature of plasma cell fate determination35, this synchrony of initiation likely resulted in a wave of differentiating plasma cell migrating from the germinal center of individuals with B cell memory.
PBMC eigenfunctions contained a signature suggesting migration of myeloid/DC lineage cells 24 hours after vaccination. Our results are quite consistent with a recent study using transcriptional analysis of whole blood after influenza vaccination5. In this study, neutrophils and monocytes were implicated as the major source of the IFN-related gene set upregulated within 15 hours post-vaccination for influenza. Interferon-related transcripts upregulated one day after Influenza vaccination in an earlier study were also correlated with anti-Influenza antibody titers6. This myeloid/DC gene signature likely triggers the initial germinal center reaction that results in plasma cell migrations later.
At the organ level in normal human subjects, we cannot sample lymph nodes and spleen. However, our analysis suggests that with higher temporal transcriptome resolution, peripheral blood does reveal evidence of prior upstream gene regulation, reflecting hidden dynamic events within lymphoid compartments. Also, upstream transcriptome analysis revealed that migration of multiple B cell populations underlie broader B cell and plasma cell surface marker phenotypes observed by flow cytometry. Further detailed experiments are needed to determine if transcriptome patterns observed are from broad changes in B cell population cell frequencies or a single dominating subpopulation. Given that our PBMC analysis reflects mixed and dynamic cell populations, any resulting gene interaction network cannot be considered strictly causal, but rather should be interpreted as a network of temporally associated transcriptome changes, likely occurring in several cell populations. Further and more detailed cell sorting work will need to be done to validate this hypothesis.
At the cellular level, division-linked B cell differentiation36,37,38 is asymmetric with key control molecules like BCL-6 bequeathed to one daughter cell39. This, along with graded IRF4 expression21, is thought to affect diversification of B cells. Our findings suggest that even with diverse B cell fates, a common transcriptomic signature can identify pre-existing vaccine-specific B cell memory, likely a marker for heterosubtypic immunity. During a pandemic, the PCgs may be useful in early prediction of both successful vaccine responses and heterosubtypic immunity.
At the individual gene level, within the PCgs we identified key plasma cell regulatory genes including those seen in other systems biology studies. TNFRSF1710 was found in both in vitro and in vivo data sets, and XBP14,40 was central to our networks. Our high-frequency sampling extends these observations with detailed time-series data capturing individual response patterns. Other genes previously associated with the post-vaccination neutralizing antibody response include KBTBD710, which did not significantly change. We also found more genes specifically associated with plasma cells than previous studies. Whether this difference is due to the statistical strength of high-frequency sampling and FPCA, or other methodological differences, will be addressed in future studies.
Finally, we would like to note that FPCA-based time-series transcriptome analysis extends beyond simple statistical analysis. The method of FPCA and the form of the results are well suited for application to machine learning in computational modeling41. In replacing instinctive or heuristic approaches with a rigorous mathematic method that can approximate behaviors of a complex biological system, FPCA lends itself to the utilization of artificial intelligence methods in models and systems biology analyses. Specifically, the resulting eigenfunctions and matrix formulations allow more rigorous characterization and modeling of complex cellular and transcriptome responses.
In conclusion, these data provide evidence for a novel transcriptomic signature associated with differentiating B cell and plasma cell migration through peripheral blood B cell pools in subjects previously vaccinated for influenza. They also suggest that with further study, temporal transcriptome patterns that accommodate individual variation may illuminate more detailed system dynamics underlying large-scale responses. The method of FPCA may also allow robust modeling and systems analysis of immune responses using artificial intelligence and machine learning methods.
This study was approved by the Research Subjects Review Board at the University of Rochester Medical Center. Informed consent was obtained from all participants. Research data were coded such that subjects could not be identified, directly or through linked identifiers, in compliance with the Department of Health and Human Services Regulations for the Protection of Human Subjects (45 CFR 46.101(b)(4)). Subject identification numbers were re-encoded for publication.
Myeloma cell lines were maintained in log-phase growth as previously described42 and used for controls in flow cytometry and ELISPOT assays: MPR-1130 (established in our laboratory), MC/CAR, Ramos cell lines (ATCC, Manassas, VA), NCI-H929 (ATCC, Manassas, VA).
Each daily blood sample for each subject was pooled from heparinized blood collection tubes and a portion removed for PBMC analysis. Lymphocytes were isolated using discontinuous Ficoll gradient centrifugation from buffy coat samples. Magnetic anti-CD235 immunoaffinity beads (Miltenyi Biotec, Auburn, CA) were used to remove RBC contaminants.
B cell samples were enriched from heparinized whole blood with RosetteSep Immunodensity separation (Stemcell Technologies, Vancouver, BC, Canada). Negative magnetic immunoaffinity bead separation columns (Miltenyi Biotec, Auburn, CA) were used with anti-CD3 (BD Bioscience, San Diego, CA) and anti-CD235 beads (Miltenyi Biotec, Auburn, CA) to further purify total B cells. Flow cytometric analysis was performed on all isolates, showing >90% purity of the isolates.
Data was collected as previously described43 with FACS Diva data acquisition software on an LSR II cytometer. Analysis of data utilized Flowjo software (Treestar, Ashland, OR). Antibodies used included anti-CD138-APC, anti-CD20-V450 (BD Bioscience, San Diego, CA), anti-CD19-PE-Texas Red, anti-CD27- PE-Cy5, (Beckman Coulter, Indianapolis, IN), anti-CD38-PE-Cy7 (eBioscience, San Diego, CA), anti-CD3-PECy5.5 (Southern Biotech, Birmingham, AL), anti-B7-FITC (Santa Cruz Biotechnology, Santa Cruz, CA), anti-CCR10-PE (Biolegend, San Diego, CA), anti-CD62L-APC-Alexa-750, and Live/Dead Aqua (Invitrogen, Grand Island, NY).
Before analysis with Flow clustering without K (FLOCK) software (https://immport.niaid.nih.gov), compensation and initial gating to remove debris and dead cells was competed using FlowJo (TreeStar). Files for each subject were then concatenated, uploaded to ImmPort and analyzed using FLOCK 1.044 to identify centroids for representative clusters across the timepoints. The following markers were used: CD19, CD20, CD138, CD3, CD27, CD38, B7, CCR10 and CD62L. Individual timepoint sample files were then analyzed against the centroids file, assigning all events to a cluster, using the Cross Sample Comparison feature of Flock. A semi supervised clustering method was used to cluster centroids to identify non-B cells (CD19−, CD3−), T cells (CD19−, CD3+), naive B cells (CD19+, CD3−, CD20+, CD27−), memory B cells (CD19+, CD3−, CD20+, CD27+), plasmablasts (CD19+, CD3−, CD20−, CD138−) and plasma cells (CD19+, CD3−, CD20−, CD138+). Percentages of each phenotype were calculated for each sample.
Hemagglutinin assays were performed as previously described45 using the 2012–2011 WHO Influenza reagent Kit for identification of influenza isolates (WHO Collaborating Center for Surveillance, Epidemiology and Control of Influenza) and the Seiken method to remove non-specific inhibitors of hemagglutination.
RNA was extracted with the Qiagen RNeasy micro kit. Concentrations were determined by UV spectrophotometry (Nanodrop) and integrity of ribosomal RNA was confirmed with the Agilent Bioanalyzer. Barcoded sequencing libraries were prepared with Illumina TruSeq RNA kits as recommended by Illumina, using 100 ng total RNA as input. All samples from an individual subject were sequenced in the same run using an Illumina Genome Analyzer IIx (22 samples per 8 lanes). CASAVA programs (Illumina, version 1.7) were used to demultiplex samples according to barcode, align sequences with the genome, and count how many sequences aligned with known variants of RefSeq genes (data are not separated according to splice variants). Sequences (reads) per gene per million total RefSeq reads were computed, with no adjustment for transcript length.
Immediately after collection, peripheral blood cell isolates were lysed in RLT buffer, passed through a Qiashredder column (Qiagen, Germantown, MD) and snap frozen in liquid nitrogen. Samples were stored at −70°C until RNA extraction, performed using the RNeasy Micro Kit with on-column DNAse (Qiagen). Quantitative RT-PCR was performed on three experimental replicates per sample in TaqMan® Array Fast Plates with TaqMan® Universal PCR Master Mix, No AmpErase® UNG on the StepOnePlus™ Real-Time PCR System (Life Technologies, Carlsbad, CA) with 10 ng cDNA into each 10 ul reaction. RNA quality was confirmed by the presence of intact rRNA with an Agilent 2100 Bioanalyzer (Santa Clara, CA). For analysis, relative quantity (RQ) values for each gene were calculated using ΔΔCT (Cq) method46. Cq values for each gene were normalized to the mean Cq from a group of six reference genes (ACTB, POLR2A, MTA2, THRAP3, USP4, CERS2) that were chosen empirically because they were unaffected by vaccination in all five subjects from whom RNAseq data were available.
Photocleavable antibody-coated paramagnetic beads were prepared as previously published23. B cell isolates were plated in ELISPOT plates coated as indicated below and incubated at 37°C for 5 hours. The plates were washed and plates for IgG and IgA detection were soaked overnight in PBS + 0.1% Tween20. Plates for IgM detection were stored at 4°C overnight in PBS before a 1 hr soak. At each time point, ELISPOT plates were developed as previously published using Alkaline Phosphate Substrate Kit III (Vector Laboratories, Burlingame CA)15. ELISPOT well membranes were removed from plates and scanned on an HP Scanjet 8300 high-resolution scanner and well images analyzed using ELIScan and EXPLORAspot custom software as previously published47.
Wells of 96-well ELISPOT and ELISA plates were coated with 3.2 mcg/ml 2010–2011 trivalent influenza vaccine for assay of TIV-specific antibody secretion for virus-specific reactivity, or 4 ug/ml goat-anti-human IgM + IgA + IgG (KPL, Gaithersburg, MD) for total IgA, IgG, or IgM. Detection reagents for separate isotype-specific assays included goat anti-human IgM, goat anti-human IgG, and goat anti-human IgA reagents that were phosphatase-labeled for ELISPOT (KPL, Gaithersburg, MD) or peroxidase-labeled for ELISA (KPL, Gaithersburg, MD). ELISA assays were developed using 1-Step Ultra TMB substrate (Pierce Biotechnology, Rockford, IL).
Significance testing of the time-series RNA Seq data was performed using Eigen-bases estimated from Functional Principal Component Analysis as published in Wu et al13. and a test statistic of the goodness of fit of the model. The multiple test correction method proposed by Benjamini and Hochberg14 was followed to control the false discovery rate (FDR).
Data sets containing RNA Seq gene expression values were uploaded into Ingenuity Pathways Analysis (IPA) (http://ingenuity.com/). Each identifier was mapped to its corresponding object in Ingenuity's Knowledge Base. Network eligible molecules were overlaid onto a global molecular network developed from information contained in Ingenuity's Knowledge Base. Networks of Network Eligible Molecules were then algorithmically generated based on their connectivity. A series of network figures with overlaid daily gene expression data was made using Upstream Analysis in IPA, Adobe Illustrator and Acrobat Professional CS5 (Adobe, San Jose, CA). Upstream analysis modulator Z-scores of 2 and above were considered probable activators (−2 for inhibitors).
We performed the standard GO functional enrichment analysis using DAVID (http://david.abcc.ncifcrf.gov/). The reported p-values were derived from the EASE score probability ‘p-Value’, and a modified Fisher's exact test that is more conservative than the standard Fisher's exact test. Benjamini-Hochberg correction was used. ‘BIOCARTA’, ‘KEGG_PATHWAY’, and ‘REACTOME_PATHWAY’ were used for pathway enrichment analysis of these clusters. We set the DAVID EASE threshold to default and ranked the enriched functions and pathways for each cluster.
Previously published in vitro microarray gene expression data23 from CpG-stimulated normal IgM- human B cells were examined for the 742 gene PCgs seen in our in vivo data. Where multiple probe sets corresponded to one gene, we reported the single probe set with the highest interquartile range.
Pearson's test for pair-wise correlation was used to test each gene and the percentage of each cellular phenotype. The Benjamini-Hochberg procedure14 was used to control FDR at 0.05 level.
For the significance test of changes in the percent of cells by phenotype and antibody levels, we used a variant of the EDGE method48 without multiple testing correction. For the day of 50% increase, we used a method similar to Zaslavsky et al49.
M.S.Z. and H.W. conceived the study. A.D.H., M.S.Z. and H.W. designed the experiments. A.D.H. performed the experiments and transcriptomic analysis and S.L.W. performed the FPCA analysis. H.W. supervised the data analysis. S.W., X.Q, M.S., H.Y., Z.L. and J.H-W. performed the statistical analysis of the data, Z.L. performed GSEA analysis, M.R. performed FLOCK analysis of flow cytometric data, S.L.W. supervised the RNA data collection and pre-analysis/processing. A.D.H., M.S.Z., S.L.W., S.W., X.Q. and H.W. wrote the manuscript.
We would like to thank Tim Mosmann, David Topham, Mark Sangster, Tim Bushnell, Danielle C. Alcéna, Shannon Hilchey, and the members of the University of Rochester Center for Biodefense Immune Modeling for spirited discussions that greatly improved the manuscript. We would like to thank study nurse Margaret McGrath as well as Yanfang Huang, Sarah Alef, Danielle Morsch, Michelle Zanche and the URMC genomics core, and the URMC Flow Core for their expert technical assistance. This work was supported by NIH grant number HHSN272201000055C (A.D.H., M.S.Z., S.L.W., J.H.-W., M.R., M.S., S.W., X.Q., Z.L., H.W.), R01 AI098112 and R01 AI069351 (M.Z., M.R., A.D.H.)