High-dimensional profiling reveals phenotypic heterogeneity and disease-specific alterations of granulocytes in COVID-19

Abstract

Since the outset of the COVID-19 pandemic, increasing evidence suggests that the innate immune responses play an important role in the disease development. A dysregulated inflammatory state has been proposed as a key driver of clinical complications in COVID-19, with a potential detrimental role of granulocytes. However, a comprehensive phenotypic description of circulating granulocytes in severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2)-infected patients is lacking. In this study, we used high-dimensional flow cytometry for granulocyte immunophenotyping in peripheral blood collected from COVID-19 patients during acute and convalescent phases. Severe COVID-19 was associated with increased levels of both mature and immature neutrophils, and decreased counts of eosinophils and basophils. Distinct immunotypes were evident in COVID-19 patients, with altered expression of several receptors involved in activation, adhesion, and migration of granulocytes (e.g., CD62L, CD11a/b, CD69, CD63, CXCR4). Paired sampling revealed recovery and phenotypic restoration of the granulocytic signature in the convalescent phase. The identified granulocyte immunotypes correlated with distinct sets of soluble inflammatory markers, supporting pathophysiologic relevance. Furthermore, clinical features, including multiorgan dysfunction and respiratory function, could be predicted using combined laboratory measurements and immunophenotyping. This study provides a comprehensive granulocyte characterization in COVID-19 and reveals specific immunotypes with potential predictive value for key clinical features associated with COVID-19.

Keywords: COVID-19; eosinophil and basophil activation; high-dimensional flow cytometry; neutrophil heterogeneity; viral immune responses.

Fig. 1.

Disease severity-dependent neutrophilia and marked decrease of eosinophils and basophils in COVID-19. (A) Experimental and analytical workflow of the study performed on samples from moderate and severe COVID-19 patients and age-matched healthy controls. (B) Dot plots (Left) and alluvial diagram (Right) describing demographics and the clinical characteristics of the patients included in the study cohort. Individual patient values are reported in SI Appendix, Fig. S1A. Severity groups (moderate, blue; severe, purple; deceased, gray) are indicated. Oxygen scale: 1, no oxygen supply; 2, oxygen < 10 L/min; 3, low flow of oxygen (noninvasive) 10 L/min to 15 L/min; 4, high flow of oxygen (noninvasive); and 5, oxygen supply via ventilator (invasive) or extracorporeal membrane oxygenation. (C) Gating strategy for the identification of granulocyte subsets and their UMAP projection is shown on Left. On Right, analog results are obtained through unsupervised subset identification based on UMAP projection of total granulocytes. (D and E) Absolute cell counts (D) and frequencies (E), based on Trucount flow cytometry analysis, among total leukocytes for neutrophils, eosinophils, and basophils in healthy controls (n = 17), moderate COVID-19 patients (n = 10), and severe COVID-19 patients (n = 16). (F) Ratios of granulocyte absolute counts over lymphocyte absolute counts in healthy controls (n = 17), moderate COVID-19 patients (n = 10), and severe COVID-19 patients (n = 16). D–F use Kruskall−Wallis test and two-stage Benjamini, Krieger, and Yekutieli test. Bars represent median. *P < 0.05; **P < 0.01; ***P < 0.001; ****P < 0.0001. BMI, body mass index; WBC, white blood cells; CRP, C-reactive protein; SOFA, sequential organ failure assessment; SSC, side scatter; FSC, forward scatter; neu, neutrophils; eos, eosinophils; bas, basophils; lymphs, lymphocytes.

Fig. 2.

Emergence of immature neutrophils in severe COVID-19. (A) Representative plots showing neutrophil identification and CD16^bright and CD16^dim subset discrimination in different study groups. Absolute cell counts (B) and frequencies (C) among total neutrophils of CD16^bright and CD16^dim cells in healthy controls (n = 17), moderate COVID-19 patients (n = 10), and severe COVID-19 patients (n = 16). (D) UMAPs on concatenated files of total neutrophils in healthy controls and moderate and severe COVID-19 patients. (E) UMAP of CD16^bright and CD16^dim cells within pooled neutrophils from all study groups. (F) Median fluorescence intensity (MFI) of selected surface markers in the UMAP described in D. (G) Hierarchical clustering of the mean z score of the MFI of the markers in CD16^bright and CD16^dim neutrophils in healthy controls (n = 17) and moderate (n = 10) and severe (n = 16) COVID-19 patients. (H) PCA based on marker expression (MFI) on CD16^bright and CD16^dim neutrophils. B and C use Kruskall−Wallis test and two-stage Benjamini, Krieger, and Yekutieli test. Bars represent median. *P < 0.05; **P < 0.01; ****P < 0.0001.

Fig. 3.

Phenograph analysis identifies activated neutrophil subsets correlating with markers of antiviral immune response in moderate COVID-19 patients. (A) Distribution of the 27 identified phenograph clusters (k-nearest neighbors, KNN = 235) overlaid on the UMAP projection. (B) Frequency of the phenograph clusters within healthy controls (n = 17) and moderate (n = 10) and severe (n = 16) COVID-19 patients; asterisks indicate statistically significant differences compared to healthy controls; # indicates statistically significant differences between moderate and severe patients. (C) Hierarchical clustering of the MFI z score of the markers in the phenograph clusters (k means = 4). Right annotation shows the log10 ratio of cluster enrichment in moderate (blue) versus severe (purple) COVID-19 patients. (D) Hierarchical clustering of Spearman correlation between phenograph clusters and soluble factors detected in all COVID-19 patients (k means = 4). Bottom annotation shows the log10 ratio of cluster enrichment in moderate (blue) versus severe (purple) COVID-19 patients. (E) Pathway analysis based on Spearman correlation values between soluble factors and phenograph cluster frequencies of group 1 from D. Only pathways displaying a normalized enrichment score of <−1 or >1 are shown. Dark blue box indicates False Discovery Rate (FDR) < 0.05. (F) Number of significant interactions between neutrophil subsets and other circulating immune cells from mild (Top) and severe (Bottom) COVID-19 patients as determined by applying the CellPhoneDB algorithm on a publicly available scRNAseq dataset (19). (G) Interaction strength (dot color) and significance (dot size) of selected ligand (pink)−receptor (blue) pairs between neutrophils and other immune subsets in mild (Top) and severe (Bottom) COVID-19 patients (19). Interactions differing between mild and severe patients are highlighted in colored boxes. Immat Neu, immature neutrophils; Mono, monocytes; Prol. Cells, proliferating cells. B uses Kruskall−Wallis test and two stage Benjamini, Krieger, and Yekutieli test. *P < 0.05; **P < 0.01; ***P < 0.001; ****P < 0.0001.

Fig. 4.

Activation of eosinophils in patients with COVID-19. (A) Overall and individual UMAPs on concatenated files of eosinophils in healthy controls and moderate and severe COVID-19 patients. (B) UMAP displaying the MFI of selected markers. (C) MFI of selected eosinophil markers in healthy controls (n = 17) and patients with moderate (n = 8) and severe (n = 14) COVID-19. Median values for each group are indicated. (D) PCA based on MFI expression of eosinophil markers. (E) CD69⁺ eosinophils from a representative healthy control and COVID-19 patients. (F) Frequencies and absolute counts of CD69⁺ eosinophils in healthy controls (n = 17) and moderate (n = 8) and severe (n = 14) COVID-19 patients. (G) Marker expression (z score) on CD69⁻ and CD69⁺ subsets. (H) MFI of selected markers in eosinophil subsets in healthy controls (n = 15) and moderate (n = 8) and severe (n = 11) COVID-19 patients. (I) Heatmaps demonstrating Spearman correlations (r < −0.4 or r > 0.4, P < 0.05) between absolute eosinophil counts/eosinophil frequency of total leukocytes/frequency of CD69⁺ eosinophils of total eosinophils and soluble factors in COVID-19 patients. (J) Pathway analysis based on Spearman values obtained from correlation between soluble factors and frequency of total eosinophils among total leukocytes (Left) or frequency of CD69⁺ eosinophils (Right). Only pathways displaying a normalized enrichment score of <−1 or >1 are shown. Dark orange box indicates FDR < 0.05. Statistical analysis by Kruskall−Wallis test and two-stage Benjamini, Krieger, and Yekutieli test are used in C, F and H, and Wilcoxon matched-pairs signed rank test is used in H. FDR adjusted P values in C, F, and H and P values < 0.05 in J are indicated. *P < 0.05; **P < 0.01; ***P < 0.001; ****P < 0.0001.

Fig. 5.

Basophils in COVID-19 patients display an activated phenotype. (A) UMAPs on concatenated files showing basophil overview and separate plots for healthy controls and moderate and severe COVID-19 patients. (B) UMAP showing the level of expression of selected markers on basophils. (C) MFI for selected markers on basophils in healthy controls (n = 17) and moderate (n = 10) and severe (n = 12) COVID-19 patients. Median values for each group are indicated. (D) PCA and biplot based on the MFI expression of basophil markers in healthy controls and moderate and severe COVID-19 patients. (E) Heatmap demonstrating Spearman correlation (r < −0.4 or r > 0.4, P < 0.05) between absolute basophil counts/basophil frequency of total leukocytes and basophil-associated soluble factors in COVID-19 patients. (F) Pathway analysis based on Spearman values obtained from correlation between soluble factors and absolute counts of basophils. Only pathways displaying a normalized enrichment score of >1 are shown. Dark blue box indicates FDR < 0.05. Significant differences between healthy controls and patient groups in C were evaluated with Kruskall−Wallis test and two-stage Benjamini, Krieger, and Yekutieli test. FDR adjusted P values in C and P values < 0.05 in F are indicated. *P < 0.05; **P < 0.01; ***P < 0.001.

Fig. 6.

Eosinophil activation and neutrophil maturation are relevant for prediction of SOFA score and respiratory function. (A) Diagonal correlation matrix including predicted clinical outcomes and the explanatory variables included in the final models. (B–G) Multivariate linear regressions for prediction of SOFA score (B and C), PaO₂/FiO₂ ratio (D and E), and peak oxygen need (F and G) were modeled based on immune traits and laboratory markers. (B, D, and F) Scatter plots of actual values versus predicted values of the linear models. (C, E, and G) Individual contribution of the variables included in the final model.

Fig. 7.

Granulocyte phenotype in convalescent patients is largely restored. (A) Schematic illustration of paired acute (Ac) and convalescent (Cv) patient samples included in the study. (B) Absolute cell counts (based on Trucount flow cytometry approach; see Materials and Methods) and (C) granulocyte frequencies over total leukocytes in paired Ac and Cv patients (n = 15). (D) PCA based on cell counts, frequencies, and normalized MFI (z score) of all measured markers. (E) Variable correlation plot. (F) Relative contribution of variables to PC1 and PC2. (G) CD16^dim neutrophil frequencies in paired Ac and Cv patient samples (n = 15). (H) Expression of selected neutrophil markers in paired Ac and Cv patients (n = 15). (I) CD69 eosinophil expression in paired Ac and Cv patients (n = 15). (J) Expression of selected eosinophil markers in paired Ac and Cv patient samples (n = 15). (K) CXCR4 basophil expression in Ac and Cv patients. (L) Expression of selected basophil markers in paired Ac and Cv patients (n = 15). Wilcoxon matched-pairs rank test was used; ns, not significant; *P < 0.05; **P < 0.01; ***P < 0.001; ****P < 0.0001. In B, C, H, J, and L, healthy controls (n = 28, median ± interquartile, IQR) ranges are shown in gray. Bars indicate statistical significance considering all sampled patients (black) or only moderate (blue) or only severe (purple) patients.