Landscape of coordinated immune responses to H1N1 challenge in humans

Abstract

Influenza is a significant cause of morbidity and mortality worldwide. Here we show changes in the abundance and activation states of more than 50 immune cell subsets in 35 individuals over 11 time points during human A/California/2009 (H1N1) virus challenge monitored using mass cytometry along with other clinical assessments. Peak change in monocyte, B cell, and T cell subset frequencies coincided with peak virus shedding, followed by marked activation of T and NK cells. Results led to the identification of CD38 as a critical regulator of plasmacytoid dendritic cell function in response to influenza virus. Machine learning using study-derived clinical parameters and single-cell data effectively classified and predicted susceptibility to infection. The coordinated immune cell dynamics defined in this study provide a framework for identifying novel correlates of protection in the evaluation of future influenza therapeutics.

Keywords: Adaptive immunity; Immunology; Infectious disease; Influenza; Innate immunity.

Figure 1. Influenza challenge model in human volunteers.

(A) Schedule of assessments performed throughout the screening, confinement, and follow-up phases of the study. Study days –1 and 1 are consecutive calendar days. Intranasal virus inoculation was performed on study day 1 immediately after blood draw. NP, nasopharyngeal. (B) Of the screened population (n = 437), the distributions of ages for those with HAI titers >1:10 (eligible for inclusion) and with HAI titers ≤1:10 (ineligible). (C) Viral titers as measured by qRT-PCR in nasopharyngeal swabs (n = 19 virus shedders). Viral titers were below the limit of detection for all participants on study day –1, and for 16 individuals throughout the study (nonshedders). (D) Stalk-specific and full-length-HA antibody seroconversion measured by ELISA. Values shown are day 29 relative to day –1. (E) Mean daily symptom score as determined by participant-reported symptom scorecard. (F–I) Baseline-normalized values for maximum daily oral temperature (F), mean pulse rate (G), systolic blood pressure (H), and diastolic blood pressure (I). Vital signs were measured up to 4 times daily. (J–M) Baseline-normalized plasma cytokine and chemokine levels, in relative fluorescence units (RFU), measured by Luminex assay for IP-10 (J), TRAIL (K), IL-10 (L), and eotaxin-2 (M). In E–M, data for virus shedders are indicated with filled squares and solid lines; data for nonshedders are indicated with open squares and dashed lines (mean ± SEM). In F–M, values plotted are normalized to baseline (average of day –1, 1). (C–M) n = 35 volunteers. Welch’s t test (B and D); Bonferroni’s adjusted P value of the time-shedding interaction term (E–M). Red arrows indicate the time point of virus inoculation throughout all figures.

Figure 2. Simplified gating strategy used to define major immune cell subsets.

Fixed whole blood was stained with a 42-marker metal-conjugated antibody panel and analyzed by mass cytometry. Simplified gating strategy used to define major immune cell subsets is shown for a single representative sample (for complete gating strategy see Supplemental Figure 5).

Figure 3. Kinetics of T cell subsets during H1N1 infection.

(A) Relative abundances of total (left) and activated (CD38⁺Ki67⁺, right) CD4⁺ T cells, naive CD4⁺ T cells (CD45RA⁺), and memory CD4⁺ T cells (CD45RA^–). (B) Relative abundances of total (left) and activated (CD38⁺Ki67⁺, right) CD8⁺ T cells, naive CD8⁺ T cells (CD45RA⁺CD27⁺), and memory CD8⁺ T cells (CD45RA^–CD27⁺). (C) Representative biaxial plots showing the extent of CD38 and Ki67 expression on CD4⁺ and CD8⁺ T cells at days 1, 4, and 8 for a representative virus shedder. In B and C, values plotted are relative abundance (percentage of CD66^– cells normalized to baseline [average of day –1, 1]) (mean ± SEM); filled squares and solid lines indicate virus shedders, and open squares and dashed lines indicate virus nonshedders. (A and B) n = 35 volunteers. Bonferroni’s adjusted P value of the time-shedding interaction term. For plots of additional gated cell populations, see Supplemental Figure 6.

Figure 4. Kinetics of B cell subsets during H1N1 infection.

(A) Relative abundance of total B cells, non–class-switched memory (NCSM) B cells (CD27⁺IgM⁺), naive B cells (CD27^–IgM⁺), and plasmablasts (CD27⁺IgM^–CD38⁺). (B) Relative abundance of p-STAT5⁺ plasmablasts (top) and representative histogram of p-STAT5 expression on plasmablasts at days 1 and 8 for a representative virus shedder (bottom). In A and B, values plotted are relative abundance (percentage of CD66^– cells normalized to baseline [average of day –1, 1]) (mean ± SEM); filled squares and solid lines indicate virus shedders, and open squares and dashed lines indicate virus nonshedders. (A and B) n = 35 volunteers. Bonferroni’s adjusted P value of the time-shedding interaction term. For plots of additional gated cell populations, see Supplemental Figure 6.

Figure 5. Kinetics of innate immune cell subsets during H1N1 infection.

(A) Relative abundances of total (first and third panels) and activated (CD38⁺Ki67⁺, second and fourth panels) populations of CD56⁺CD16^– NK cells and CD56^loCD16⁺ NK cells. (B) Relative abundances of pDCs, CD1c⁺ myeloid DCs (mDCs), and BDCA3⁺ mDCs. (C) Relative abundances of classical monocytes (cMCs; CD14⁺CD16^–), intermediate monocytes (intMCs; CD14⁺CD16⁺), and nonclassical monocytes (ncMCs; CD14^–CD16⁺). (D) Relative abundances of CD66⁺ cells and CD123⁺HLA-DR^– cells. (E and F) Biaxial plots of CD38 and Ki67 expression on CD56^loCD16⁺ NK cells (E) and CD14 and CD16 expression on monocytes (F) on days 1–8 and 29. Plots for 1 representative virus shedder are shown. In B–E, values plotted are relative abundance (percentage of CD66^– cells normalized to baseline [average of day –1, 1]) (mean ± SEM); data for virus shedders are indicated by filled squares and solid lines, and data for nonshedders are indicated with open squares and dashed lines. (A–F) n = 35 volunteers. Bonferroni’s adjusted P value of the time-shedding interaction term. For plots of additional gated cell populations, see Supplemental Figure 6.

Figure 6. Time-based Scaffold analysis of volunteer H1N1 challenge.

Pooled cells from all volunteers and all time points (381 mass cytometry samples) were binned into 1 of 200 nonoverlapping clusters based on the similarity of their expression across all surface markers. Clusters are displayed as nodes of a force-directed Scaffold graph with the length of edges inversely proportional to the similarity of clusters to each other and to landmark nodes defined by manual gating (Supplemental Figure 5). The colored wedges in cluster nodes indicate the degree of increase (red) or decrease (blue) in the mean relative abundance of a given cell cluster for virus shedders on study days 2–8. n = 35 volunteers.

Figure 7. Kinetics of Scaffold cell clusters during H1N1 infection.

(A) Nodes colored black on the Scaffold graph (see Figure 6) indicate clusters with significant differential abundance between shedders and nonshedders (P < 0.05, Bonferroni’s adjusted P value of the time-shedding interaction term). The locations of the 7 Scaffold clusters selected for further analyses are indicated by Roman numerals. (B) Line plots showing abundances (percentage of total clustered cells) (mean ± SEM) normalized to baseline (average of day –1, 1) in Scaffold clusters indicated in A. Data for virus shedders are indicated with filled squares and solid lines, and data for nonshedders are indicated with open squares and dashed lines. (C) Box plots showing expression levels of indicated cell surface markers in Scaffold clusters VI and VII relative to manually gated myeloid and plasmacytoid DC landmarks. Note differential expression of CD38. (A–C) n = 35 volunteers. Bonferroni’s adjusted P value of the time-shedding interaction term. For additional box plots showing all markers for Scaffold clusters I–VII, see Supplemental Figure 7.

Figure 8. Analysis of the role of CD38 in pDC activation.

(A) Median CD38 values on manually gated pDCs analyzed by mass cytometry. Baseline-normalized (average of day –1,1) values (mean ± SEM) are plotted for virus shedders (filled squares and solid lines) and nonshedders (open squares and dashed lines). (B) Histograms of CD38 expression on pDCs at days 1, 3, 5, and 8 for a representative virus shedder. (C and D) PBMCs from healthy donors were cultured for 24 hours in the presence or absence of the indicated TLR agonist and analyzed by flow cytometry. (C) Representative histograms for CD38 expression on pDCs (see pDC gating strategy in Supplemental Figure 9). (D) CD38 MFI values from experiments with cells pooled from 4–11 donors in 2–7 independent experiments. (E–K) PBMCs were cultured for 24 hours in the presence (H1N1, blue) or absence (No stim., orange) of heat-inactivated H1N1 virus and, where indicated, in the presence of anti-CD38 blocking antibody or isotype control (Isotype). Flow cytometry data from PBMCs from 3 donors in 2 independent experiments are shown. (E) Representative histograms of CD38 expression on pDCs. (F) Quantification of CD38 MFI on pDCs. (G) Representative biaxial plot showing frequency of TNF-α⁺ pDCs. (H) Percentage of TNF-α⁺ pDCs. (I) Representative biaxial plots showing frequency of IFN-α⁺ pDCs. (J) Percentage of IFN-α⁺ pDCs. (K) Soluble IFN-α measured in supernatants of PBMC cultures. Bar plots in D, F, H, J, and K show mean ± SEM. Welch’s t test (D, F, and H); Wilcoxon’s signed-rank test (J and K).

Figure 9. Machine learning–based classification and prediction of H1N1 virus shedding status.

(A) Schematic of models used to classify virus shedding after challenge. Four models were generated from 4 data sets: cytokines/chemokines; CBC; cell abundance; and cell activation and proliferation (Supplemental Table 8). For each data set, the model was iteratively trained and tested on 50% held-out data 100 times. The performance for each individual model was the median value from 100 iterations. For the stacked model, individual models were combined. (B) Heatmap of P values for classification of virus shedding status on days 2–8. (C) Correlation network displays the relationship between features used in the combined model. Node sizes are proportional to the univariate correlation between a given feature and virus shedding. (D) The P values for each model on the day with the best overall classification power (day 6). (E) Receiver operating characteristic (ROC) curve evaluating the performance of the combined classification model on day 6. (F) Schematic of model used to predict virus shedding before challenge. The 4 data sets from A were used to generate a predictive model that was trained and tested iteratively on 50% held-out data 100 times. In contrast to A, only the top features associated with classification of virus shedding on day 6 from each data set were used to train models on day 1 data. Models were subsequently evaluated on the excluded day 1 test set. The performance for each individual model was the median value from 100 iterations. For the final model, individual models were combined. (G) The P values for each model on day 1. (H) ROC curve evaluating the performance of the combined predictive model on data collected on day 1. Wilcoxon’s signed-rank test (A–H). n = 35 volunteers.