Mitochondrial induced T cell apoptosis and aberrant myeloid metabolic programs define distinct immune cell subsets during acute and recovered SARS-CoV-2 infection

Main Text

SARS-CoV-2 is a coronavirus responsible for the COVID-19 pandemic, resulting in over 20 million cases worldwide. While the vast majority of infected patients experience a self-limiting viral syndrome, others develop severe disease leading to pneumonia and acute respiratory distress syndrome (ARDS), which has accounted for over 900,000 deaths globally (1). At this time, it is unclear why some patients readily resolve infection while others develop severe symptoms. Specifically, it remains to be determined if severe disease is associated with a failure to generate protective immunity, overly robust dysfunctional immune responses, or a combination of both.

It has become increasingly clear that metabolic reprogramming is not just a consequence of immune activation, but rather plays a critical role in facilitating immune cell differentiation and function (2, 3). It is now understood that discrete cellular subsets employ distinct metabolic programs. For example, effector T cells are characterized by increased expression of molecules necessary to support glycolysis, while memory T cells upregulate expression of molecules involved in oxidative phosphorylation and fatty acid oxidation (4). Exhausted T cells are characterized not just by the upregulation of inhibitory molecules such as programmed cell death protein 1 (PD-1) and loss of polyfunctionality, but also by mTOR signaling in the absence of productive glycolytic function and anabolic processes (5). Therefore, combining immunologic markers with metabolic markers has the potential to transcend traditional immune cell phenotyping and provide novel insights into distinct functional subsets. To achieve this, we created unique high-dimensional (HD) flow cytometry panels that combine lymphocyte and myeloid markers with those of metabolic function and used them to examine peripheral blood mononuclear cells (PBMC) from acutely infected COVID-19 patients (COVID-A) (Table S1, Fig. S1).

HD flow cytometry was performed on thawed PBMC from an IRB approved biorepository from patients admitted to the Johns Hopkins Hospital (Table 1). We initially focused on T cells within the PBMC given their importance in viral control. While we did not observe significant differences between the percentages of CD3⁺, CD4⁺, or CD8⁺ T cells between the PBMC from COVID-A patients and healthy controls, we observed an increase in the CD4:CD8 ratio in COVID-A patients, as previously reported (Fig. S2) (6, 7). Otherwise, conventional T cell markers simply revealed an increase in central memory T cells (Tcm) in the CD4⁺ T cell population of the COVID-A patients when compared to healthy controls.

In stark contrast, unbiased analysis of T cells employing the combined immune and metabolic markers revealed the robust presence of a unique subset of both CD8⁺ and CD4⁺ T cells in the PBMC from COVIDA patients (Fig. 1A). Having identified this unique population of T cells, we next sought to determine their precise phenotype. When compared to healthy controls, we observed no differences in expression of the classical markers of T cell subsets CD45RA, CCR7, or KLRG1 (Fig. S2F) or the activation markers CD69, Ki67, PD-1, or HLA-DR (Fig. S2G). Likewise, levels of metabolic enzymes involved in glycolysis and fatty acid oxidation did not differ between the two groups (Fig. S2H). In contrast, the unique population of T cells found in the COVID-A patients was characterized by robust upregulation of Voltage Dependent Anion Channel (VDAC) and the epigenetic mark H3K27me3 (Fig. 1A-B). H3K27me3 is regulated in part by a-ketoglutarate-mediated jumonji demethylases (8), while VDAC is a mitochondrial membrane protein involved in metabolite transport and has been associated with promoting cell death and lupus-like autoimmunity (9). Upon typical T cell activation, as seen upon TCR stimulation of healthy PBMCs in vitro, VDAC and H3K27me3 expression increase along with the glucose transporter Glut1 and Hexokinase II (HKII), both of which support activation-induced glycolysis (Fig. 1C). However, the unique population of H3K27me3^hiVDAC^hi T cells in the PBMC of COVID-A patients do not appear to resemble conventional recently activated T cells because they express relatively low levels of Glut1 and HKII (Fig. 1D). Furthermore, the H3K27me3^hiVDAC^hi T cells encompass both CD45RA/CCR7 positive and negative cells (Fig. 1D) and did not demonstrate enrichment for specific TCR clones (Fig. 1E). Expanded analysis of 55 PBMC samples from 38 COVID-A patients (including sequential samples) revealed that a significant number of patients possess high proportions of T cells with this phenotype (Fig. 1F). Interestingly, while not present in all of the patients queried, every patient 70 years of age and older possessed increased frequencies of the H3K27me3^hiVDAC^hi T cells in their PBMC (Fig. 1G). Overall, these data reveal a novel, distinct population of H3K27me3^hiVDAC^hi T cells in the PBMC of COVID-A patients. The discordance in activation and metabolic programs suggests that these T cells are dysfunctional.

• Download figure
• Open in new tab

Figure 1: Identification of novel metabolically distinct T cells in COVID-19 patients

A. Concatenated flow cytometry data depicted as UMAP projection of CD3⁺ T cells from healthy control (HC, grey) and hospitalized COVID-19 patients (COVID19, blue). The two markers discovered to drive segregation of the COVID and HC cluster, H3K27Me3 and VDAC, are depicted as histogram overlays and MFI heatmap overlays on UMAP projection. B. Representative flow plots of H3K27Me3⁺VDAC⁺ T cells. C.Healthy PBMCs were stimulated for 24 hours with anti-CD3 and anti-CD28 and evaluated for metabolic enzymes using flow cytometry. Representative histograms of unstimulated (grey) and stimulated (purple) cultures. D. Histograms comparing H3K27Me3⁺VDAC⁺ T cells (pink) and remaining T cells (blue) from concatenated pooled COVID-19 donors for indicated proteins. E. TCR Simpson Clonality from sorted H3K27Me3⁺VDAC⁺ T cells (pink) and remaining T cells (blue). F. Frequency of H3K27Me3⁺VDAC⁺ as percent of CD4 or CD8 T cells. Each dot represents one individual, significance tested using unpaired Mann-Whitney test. G. Frequency of H3K27Me3⁺VDAC⁺ as percent of CD4 (black) or CD8 (grey) stratified by age of COVID-19 patients. H. UMAP projection of pooled donors with active infection, color coded by disease. I. UMAP projection of active and recovered COVID19 compared to healthy controls. J. UMAP projection of influenza and active COVID-19 compared to healthy controls. K. Frequency of H3K27Me3⁺VDAC⁺ CD4 and CD8 T cells as percent of total live cells. L. Representative gating of H3K27Me3⁺VDAC⁺ CD4 T cells from multiple disease states. M. Hierarchical clustering of H3K27Me3⁺VDAC⁺ CD4 T cells based on expression (MFI values) of indicated proteins. Comparison of severe hospitalized COVID-19 infection (dark blue), recovered COVID-19 (light blue) and hospitalized influenza (red). N. Normalized MFI of GLUT1, TOMM20 and KLRG1 in all patients and representative histogram overlays of MFI. Each dot represents one individual, significance tested using unpaired Mann-Whitney test (F), two-way ANOVA (G) or unpaired Kruskal-Wallis test compared to healthy control (K) or every combination (N). *p < 0.05, **p < 0.01, ***p < 0.001, and ****p < 0.0001

We then tested whether the unique population of H3K27me3^hiVDAC^hi T cells in COVID-A patients were present in the PBMC of recovered COVID-19 patients (>28 days from diagnosis, COVID-R) and patients with other viral infections. Thus, we analyzed PBMC from patients with acute and chronic hepatitis C infection (10), hospitalized with influenza infection (11), and those who had recovered from COVID-19 (COVID-R)) (Table S2) (12). There appeared to be only subtle differences between these groups in CD4⁺ and CD8⁺ T cell subsets as defined by traditional markers (Fig. S2). Yet, when compared to the T cells in the PBMC from other groups, the COVID-A patients segregated as shown by UMAP projections, indicating a distinct subset of T cells unique to COVID-A patients not accounted for by traditional markers (Fig. 1H-J). While we did not observe the H3K27me3^hiVDAC^hi T cells in the PBMC of patients with either acute or chronic hepatitis C infection, H3K27me3^hiVDAC^hi T cells were present in the PBMC of some of the hospitalized influenza patients and of some COVID-R patients (Fig. 1K-L). For eight of the COVID-A patients, we were also able to obtain PBMC during the resolution phase, approximately 90 days after the onset of infection. Interestingly, the three patients with the highest levels of the H3K27me3^hiVDAC^hi T cells while hospitalized demonstrated markedly diminished percentages of these cells at day 90 (Fig. S3A-B).

Unbiased hierarchical analysis based on protein expression of H3K27me3^hiVDAC^hi CD4⁺ T cells showed that these cells robustly clustered based on disease type, indicating qualitative differences in H3K27me3^hiVDAC^hiT cells across disease status (Fig. 1M). That is, the H3K27me3^hiVDAC^hiT cells in the PBMC of the influenza and COVID-R patients, were distinct from the H3K27me3^hiVDAC^hiT cells in the PBMC of the COVID-A patients. For example, the H3K27me3^hiVDAC^hi T cells from the COVID-A patients demonstrated significantly decreased expression of Glut1, a metabolic marker associated with T cell effector function (Fig. 1N). However, they had significantly increased expression of the mitochondrial protein TOMM20 and KLRG1, a marker associated with T cell senesence and age-related functional defects (13), when compared to the H3K27me3^hiVDAC^hi T cells from both the Influenza patients and the COVID-R patients (Fig. 1N). The H3K27me3^hiVDAC^hi T cells present in both the influenza patients, as well as the COVID-R patients, displayed increased markers of glycolysis and effector function, and were metabolically distinct from the COVID-A patients. Thus, despite similar levels of H3K27Me3 and VDAC expression, the metabolic programming of these cells is distinct between acute and recovered COVID-19 patients and from the hospitalized influenza patients, all of whom recovered.

The percentage of H3K27me3^hiVDAC^hi T cells in the PBMC of many of the COVID-A patients exceeded 50% of total T cells (Fig. 1F). In order to elucidate the functional significance of these cells, single-cell RNA sequencing was performed on six COVID-A patients, which all had above 50% of the T cells of interest as determined by flow cytometry, and three healthy controls. Evaluating gene programs enriched in T cells derived from COVID-A patients revealed an upregulation of apoptotic signatures and cell death, as has been reported previously (Fig. 2A, red bars) (14). Furthermore, these data demonstrated evidence of mitochondrial dysfunction, including downregulation of several programs associated with mitochondria function, mitochondrion organization, mitochondrial respiratory chain complex assembly, oxidative phosphorylation, and electron transport coupled proton transport (Fig. 2A, blue bars). To investigate whether elevated VDAC and TOMM20 indicate defective mitochondrial function in these T cells, we performed electron microscopy on PBMC from additional COVID-A patients and healthy controls to examine the lymphocyte mitochondria. This analysis revealed markedly dysmorphic, irregularly shaped mitochondria with incomplete cristae in lymphocytes from COVID-A patients compared to lymphocytes from healthy controls, consistent with dysregulated mitochondria function (Fig. 2B). In addition, lymphocytes from COVID-A patients showed prominent degenerative changes including large cytoplasmic vacuoles. These observations were concurrent with cells in COVID-A PBMC displaying morphological characteristics of apoptosis. We then performed confocal microscopy of PBMC stained with anti-CD3 and MitoTracker Deep Red Dye, which stains mitochondria. CD3⁺ T cells from COVID-A patients demonstrated less distinct mitochondrial staining versus healthy controls (Fig. 2C). Having observed decreased mitochondrial programs by RNA-seq, presence of apoptotic cell morphologies by EM and confocal Deep Red Dye staining, we hypothesized that cytochrome c was being released from the mitochondria leading to apoptosis (15). We therefore performed immunofluorescence analysis of endogenous cytochrome c, and found it to be present in the cytoplasm of CD3⁺ T cells from COVID-A patients, whereas in healthy controls, cytochrome c was localized to the mitochondria (Fig. 2D). Thus, both electron and confocal microscopy support the single-cell RNA sequencing data demonstrating that the unique population of T cells found in the COVID-A patients display mitochondrial dysfunction consistent with cytochrome c release into the cytoplasm and subsequent induction of apoptosis.

• Download figure
• Open in new tab

Figure 2: Dysregulated T cell mitochondria associated with apoptosis can be rescued by targeting VDAC

A. Single-cell RNA sequencing analysis of 6 COVID-19 subjects and 3 HC were evaluated for CD3+ T cells. Genes distinguishing T cells from COVID-19 patients compared to HC were evaluated for statistical over representation using GO biological processes as gene sets and categorized into higher level annotation using ReviGO. Displayed in the enrichment score for each gene set and color corresponds to programs in upregulated genes (red) and downregulated genes (blue). B. Representative electron microscopy images of PBMCs from a COVID-19 patient and a healthy control. C. Representative confocal images of PBMCs from a COVID-19 patient and healthy control with mitochondria labeled using MitoTracker Deep Red (pink), CD3⁺ T cells labeled (green) and nuclei labed with DAPI (blue). D. Representative fluorescence images of PBMCs from 3 COVID-19 subjects and one healthy control (left) immunostained for cytochrome c (green) and CD3 (red), and nuclei labeled with DAPI (blue). Plot profiles of intracellular cytochrome c fluorescence intensity distribution (right). E. PBMCs from COVID-19 patient or healthy control were cultured for 48 hours in media, rapamycin (100nM), VBIT-4 (300nM), or zVAD (60nM). T cell survival was calculated as the percent CD4 or CD8 T cells remaining from initial plating. Significance tested using two-way ANOVA with each drug compared to the media control, n = 9. F. T cells were stimulated with anti-CD3/CD28 (purple) for 48 hours and surviving T cells from COVID-19 patients are able to respond by upregulating HLA-DR, CD69, CD25 and GLUT1 to the same extent as healthy controls compared to unstimulated controls (grey). *p < 0.05, **p < 0.01, ***p < 0.001, and ****p < 0.0001

Given these findings and the fact that VDAC facilitates caspase-mediated cell death (16), we hypothesized that the high expression of VDAC was directly linked to increased susceptibility to cell death in these T cells. To test this hypothesis, we cultured PBMC from COVID-A patients and healthy controls for 48 hours in vitro in the presence of media alone, the mTOR inhibitor rapamycin, the VDAC oligomerization inhibitor VBIT-4 (9, 17), and the global caspase inhibitor ZVAD (Fig. 2E). We observed decreased survival in media alone of the T cells from the PBMC of COVID-A patients when compared to healthy controls. Interestingly, survival of the COVID-A T cells was rescued with both the VDAC oligomerization inhibitor and the pan-caspase inhibitor, but not the mTOR inhibitor. This suggests that VDAC is promoting cell death in these T cells. Further, neither VBIT-4 nor ZVAD enhanced survival of the healthy T cells, demonstrating that this mechanism of cell death is not active in the cultured healthy T cells. In light of the propensity of the T cells from the COVID-A patients to undergo apoptosis, we wanted to test that our results did not reflect differences in cell survival after freezing. When we compared fresh and frozen PBMC from four COVID-A patients we observed an increased frequency of the H3K27me3^hiVDAC^hi T cells in the freshly stained PBMC (Fig. S3C-D) compared to those frozen and thawed. Thus, our studies are potentially underestimating the frequency of these unique T cells in the COVID-A patients. For 12 COVID-A patients with multiple sequential samples, we did observe a correlation between the presence of the H3K27me3⁺VDAC⁺ T cells and subsequent development of lymphopenia (preliminary data, Fig. S3E) supporting the hypothesis that mitochondrial-induced apoptosis might be contributing to the lymphopenia observed in COVID-A patients (6, 18). Interestingly, the COVID-A T cells, which were rescued by the ZVAD and VBIT-4 treatment, respond similarly to anti-CD3 + anti-CD28 stimulation to the healthy T cells (Fig. 2F). Specifically, stimulation of the COVID-A T cells rescued from cell death resulted in robust, productive activation.

In contrast to T lymphocytes, we did not observe marked differences in the frequencies of naïve B cells between any viral infections and the healthy controls (Fig. S4). Consistent with recent reports, we did observe an increase in the peripheral antibody-secreting cells (ASC) in the COVID-A patients (18). All of the viral infections were associated with a significant increase in the memory B cells versus healthy controls, while only the COVID-A and hepatitis C patient groups demonstrated an increase in the activated memory B cells. The COVID-A and influenza patient groups both had a higher percent of atypical memory B cells compared to PBMC from the hepatitis C, and the COVID-R patient groups. Global UMAP projection demonstrated only subtle phenotypic differences among the different groups. When we examined the NK cells, we found no significant differences in frequency when comparing the COVID-A patients, the COVID-R patients and the influenza patients (Fig. S5). Conversely, HD flow analysis revealed a population of CD56+ NK cells in the COVID-A patients not present in the influenza patients or in the healthy controls. This population was defined by the upregulation of classical activation markers CD69, Ki67, and CD49a and by increased expression of the metabolic markers TOMM20, CPT1a, and HKII. While the precise significance of these cells is unclear, their generation can potentially provide important clues into the dysregulated systemic inflammation characteristic of acutely ill COVID-A patients. When we compared NK cells from COVID-A patients with COVID-R patients, we observed differences in CD56+ and CD56bright cells driven mainly by differential expression of markers of activation.

Next, we examined myeloid cells in the PBMC of the COVID-A patients employing our HD immune-metabolic panel. There was not a difference in the frequency of total myeloid cells between all the groups examined (Fig. S6). We did observe a significant decrease in both myeloid dendritic cells (mDC) and plasmacytoid dendritic cells (pDC) in the PBMC from the COVID-A patients when compared to PBMC from the COVID-R and the influenza patient groups (Fig. S6C). Consistent with prior reports, acute COVID-19 was associated with a decreased percentage of pDC in PBMC (18, 19).

Visualization of the data by UMAP projection revealed two distinct myeloid populations, which once again became apparent only by including the metabolic markers in our panel (Fig. 3A-B). We first identified CD15⁺ granulocytic cells in COVID-A patients, which were entirely absent from the healthy controls (Fig. 3A-B). Evaluating the granulocytic cells further revealed a combination of low-density neutrophils and polymorphonuclear (PMN)-myeloid derived suppressor cells (MDSC) (Fig. 3C). Interestingly, a specific subset of the PMN-MDSC also expressed elevated levels of the rate limiting enzyme of glycolysis HKII. Increased neutrophil counts have been previously observed in COVID-19 patients (18, 19). Low-density neutrophils can be found in the buffy coat during inflammation and represent both immature neutrophils and activated/degranulated neutrophils (20). Low-density neutrophils were present in a subset of influenza patients, but only significantly enriched in the COVID-A patients (Fig. 3D). Strikingly, while the immunosuppressive PMN-MDSC (as defined by conventional markers) are detectable in both chronic hepatitis C and COVID-A patients, the HKII⁺ PMN-MDSC were unique to COVID-A patients (Fig. 3E). Thus, our approach revealed a novel population of HKII⁺ PMNMDSC which were found exclusively in the COVID-A patients.

• Download figure
• Open in new tab

Figure 3: Metabolically distinct immunosuppressive myeloid derived suppressor cells in PBMC of COVID-19 patients

A. Concatenated flow cytometry data depicted as UMAP projection of CD3⁻CD19⁻CD56⁻myeloid cells from healthy control (HC, grey) and acute COVID-19 patients (COVID, blue). B. UMAP projection of MFI heatmap overlays of indicated proteins. C. Representative gating of CD15⁺ PMN-MDSC. D-E Frequency of indicated cell subset as percent of total live cells. F. Representative gating of CPT1a⁺VDAC⁺ myeloid cells (gated on CD3⁻CD19⁻CD56⁻and CD33⁺) and subset of CPT1a⁺VDAC⁺ cells based on HLA-DR expression. GH. Frequency of HLA-DR⁻ CPT1a⁺VDAC⁺ myeloid cells and HLA-DR⁺ CPT1a⁺VDAC⁺ myeloid cells as percent of total live cells. Each dot represents one individual, significance tested using unpaired Kruskal-Wallis test compared to healthy control. *p < 0.05, **p < 0.01, ***p < 0.001, and ****p < 0.0001

When interrogating the monocytic cell populations in the PBMC, again the metabolic markers helped distinguish a unique population of cells in the COVID-A patients when compared to healthy controls. (Fig. 3F). For comparison, we did observe differences between the myeloid cells of the influenza patients and of the healthy controls; however, these distinctions were driven primarily by classical activation markers such as CD86 and CD40 (Fig. S6D). In contrast, the monocytic cells found in COVID-A patients expressed high levels of carnitine palmitoyltransferase 1a (CPT1a), an enzyme found within the mitochondrial membrane that is essential for fatty acid oxidation. In addition they expressed high levels of VDAC1 (Fig. 3F). CPT1a has been associated with both NLRP3 inflammasome activation and ROS production (21, 22). A potential role for inflammasome activation as contributing to the pathogenesis of SARS-CoV-2 infection has been noted by others (23). A large portion of the CPT1a⁺VDAC⁺ myeloid cells were CD14⁺HLA-DR⁻ and thus classified as monocytic MDSC (M-MDSC) (24). Furthermore, while these cells were significantly upregulated in the COVID-A patients, they were only minimally present in the PBMC of the COVID-R patients (Fig. 3G). A decrease in HLA-DR expression on myeloid cells during COVID-19 infection has been previously noted by other groups as well (23, 25, 26). Interestingly, the HLA-DR⁺ component, composed of both CD14⁺and CD16⁺ monocytes, was elevated during active infection and remained high during recovery (Fig. 3H). HLA-DR⁺ monocytes are capable of producing robust amounts of cytokines and also providing antigen presentation and T cell costimulation. The differential expression patterns of these monocytic populations between the infected and recovered patients may represent a shift from immunosuppressive (DR-, MDSC) to a productive immune response (DR+, stimulatory monocytes).

The dichotomy based on HLA-DR expression of these unique cells prompted us to examine the relationship between the HLA-DR⁺ and HLA-DR⁻ CPT1a⁺VDAC⁺ myeloid cells and disease severity. Not only do the CPT1a⁺VDAC⁺ HLA-DR⁻ cells represent less than 0.5% of total PBMCS of COVID-R patients, but the percentage of these cells is significantly higher in COVID-A patients requiring mechanical ventilation (on average 3.5% of PBMCs) than in COVID-A patients with less severe disease (0.7% of PBMCS) (Fig. 4A). Thus, not only can the CPT1a⁺VDAC⁺HLA-DR⁻ cells be employed as a potential biomarker of disease, but comparisons between the DR⁺ and DR⁻ subsets of these cells provides potentially important insight into the pathogenesis of disease.

• Download figure
• Open in new tab

Figure 4: Presence of immune cells with distinct metabolic profiles predicts disease severity

A. Frequency of CPT1a⁺VDAC⁺ DR^+/− cells in COVID-19 patients were stratified by disease severity (Severe = deceased or on ventilator, Mild = hospitalized with low or high flow oxygen, Recovered = post resolution of COVID19). Each dot represents one individual, significance tested using two-way ANOVA comparing either HLA-DR⁺ vs HLA-DR⁻ within each category or HLA-DR⁺ or HLA-DR⁻ across each category. B. UMAP projection of scRNA seq of myeloid cells from three COVID patients with detectable CPT1a⁺VDAC⁺ myeloid cells by flow cytometry colored by identified clusters 0–3. C. Expression of indicated gene within cluster 0–3. Each dot represents a single cell. D. Genes identifying cluster 1 (CPT1a⁺VDAC⁺ HLA-DR^dim) and cluster 3 (CPT1a⁺VDAC⁺ HLA-DR^high) were evaluated for statistical over representation using GO biological processes as gene sets and categorized into higher level annotation using ReviGO. Heatmap color corresponds to the enrichment score in upregulated genes (red) and downregulated genes (blue), x indicates a non-significant enrichment. E. Expression of indicated gene within cluster 0–3. Each dot represents a single cell. F-I. Feature importance of distinguishing COVID-19 patients and healthy controls, COVID-19 patients and Flu patients, severe COVID-19 patients and recovered COVID-19 patients, or severe COVID-19 patients and mild COVID-19 patients as indicated. Each feature is depicted as a frequency of total live cells, unless otherwise indicated in the case of proportion of cell subset with a specific metabolic phenotype.

To better understand the function of these CPT1a⁺VDAC⁺ HLA-DR⁻ cells, we assessed single-cell RNA sequencing data from the PBMC of three COVID-A patients that had very high levels of the CPT1a⁺VDAC⁺ myeloid cells, as assessed by flow cytometry (Fig. 4B). The analysis revealed four distinct clusters (Fig. 4B). Clusters 1 and 3 both had elevated levels of VDAC and CPT1a, but cluster 1 expressed relatively lower levels of HLA-DR compared to cluster 3 (Fig. 4C). Cluster 3 exhibited gene programs associated with a productive immune response, such as antigen processing/presentation and a type I IFN response (Fig. 4D). In contrast, cluster 1, which demonstrated decreased HLA-DR expression and contained the cells that were more abundant in severe disease, expressed gene programs associated with ROS production, exocytosis and targeting proteins to the membrane surface (Fig. 4C). Specifically, the HLADR⁻ cluster expressed high levels of genes for secreted alarmins such as S100A9 and S100A8 and chemokines such as CXCL2 and CXCL3 (Fig.4E) (27). Thus, not only do the HLA-DR⁻ and HLA-DR⁺ CPT1a⁺VDAC⁺ myeloid cells distinguish disease severity and recovery, but these cells express markedly different gene expression profiles that have the potential to provide important insight into disease pathogenesis and productive immune responses.

We pursued the potential role of these novel immune cells in predicting disease by treating the percentage of each cell population as features. We trained random forest models for classifying patients into different groups (e.g., COVID-A vs. Healthy controls, Severe COVID-A vs. COVID-R, Severe COVID-A vs. influenza, Severe COVID-A vs. Mild COVID-A). Consistent with our previous observations, cell populations such as H3K27Me3⁺VDAC⁺CD4⁺, H3K27Me3⁺VDAC⁺CD8⁺, and PMN-MDSC cells were among the most important features in distinguishing COVID-A patients and healthy controls (Fig. 4F). Comparing severe COVID-A patients to influenza patients also revealed cell populations, such as PMN-MDSC, mDC, and MMDSC cells, that distinguish between these two diseases (Fig. 4G). When severe COVID-A patients were compared to recovered (Fig. 4H) or mild COVID-A patients (Fig. 4I), cell populations such as pDC, M-MDSC and HKII⁺PMN-MDSC were found to be among the top 10 distinguishing features in both comparisons (Fig. 4H-I). Additionally, cell populations such as VDAC⁺CPT1a⁺ myeloid cells also had distinguishing power when comparing severe and mild COVID-A patients (Fig. 4I), whereas cell populations such as mDC were among the top features to distinguish severe and COVID-R.

In terms of prediction accuracy, the receiver operating characteristic (ROC) analysis based on leave-oneout cross-validation (LOOCV) of the random forests indicates that COVID-A patients are highly distinguishable from the healthy controls, and the severe COVID-19 patients are also highly distinguishable from the influenza or COVID-R patients, whereas distinguishing between severe and mild COVID-19 was more challenging (Fig. S7A). Because prior publications have shown the association between COVID-19 severity and age, sex, and BMI, we further added these variables to our COVID-19 severity analysis (6). The feature importance analysis confirmed in our dataset that sex and BMI are among the most important features in predicting COVID-19 severity, whereas age was less important in the presence of our cell population features (Fig. S7B). The percentage of VDAC⁺CPT1a⁺ myeloid cells, pDC, and H3K27Me3⁺VDAC⁺CD4⁺ cells were more important than sex and BMI in predicting disease severity. Based on these observations, we trained two additional prediction models using the top-fiveranked features (i.e., percentage of VDAC⁺CPT1a⁺, pDC, and H3K27Me3⁺VDAC⁺CD4⁺cells, sex, and BMI) and the basic clinical information (i.e., age, sex, and BMI), and compared their performance in predicting COVID-19 severity (Fig. S7C). The model based on the top-five-ranked features outperformed the model based on the basic information, indicating that the percentage of VDAC⁺CPT1a⁺ myeloid cells, pDC, and H3K27Me3⁺VDAC⁺CD4⁺improved prediction of COVID-19 severity.

Overall, our study has identified previously undescribed populations of immune cells that are distinctly upregulated in the PBMC of acutely infected COVID-19 patients. Future focused investigation of these unique subsets provides a platform to dissect and identify the pathogenic inflammatory modulators of disease as well as discern the mechanistic basis for the development of severe lung disease that leads to death. Indeed, the H3K27me3^hiVDAC^hi T cells demonstrate mitochondrial dysfunction leading to cytoplasmic cytochrome c release and caspase dependent apoptosis that can be inhibited by blocking VDAC oligomerization (Fig. 2). This dysregulated T cell immunity could contribute to a lack of or waning protective immunity or impair the functionality of pre-exisiting cross-reactive T cell immunity (28, 29). Further, this mitochondrial mediated cell death might contribute to the lymphopenia that is observed in COVID19 patients. Alternatively, it is possible that increased cell death of T cells is driving some of the dysregulated inflammation and autoimmune-like features characteristic of COVID-19. Finally, our data suggests that therapeutically targeting mitochondrial (metabolic) dysfunction might represent a successful strategy for abrogating disease.

Previous studies have demonstrated increased inflammation and activation of myeloid cells in COVID-19 patients (18, 23, 26). However, it still remains unclear which of these responses are protective and which are pathogenic. Our identification of unique myeloid populations in the COVID-A patients shed light on the role of these cell populations. While PMN-MDSC were observed in multiple viral infections, the HKII⁺ PMN-MDSC were exclusively upregulated in the COVID-A patients, suggesting that this population of cells is either contributing to pathogenesis or the consequence of dysregulated inflammatory responses (Fig. 3). CPT1a⁺VDAC⁺DR⁺ monocytic cells were found in both the acutely and recovered patients. DR⁺ monocytes have previously been shown to be associated with productive immune responses in other respiratory viral infections (30) and single-cell RNA sequence analysis revealed increased antigen processing/presentation and type I IFN responses in such cells. In contrast, the CPT1a⁺VDAC⁺DR⁻M-MDSC were found exclusively in the COVID-A patients and further, the percent of these cells was positively correlated with severity of disease (Fig. 3–4). As such, the presence of these cells seems to be indicative of dysregulated inflammation. Thus, these cells represent potentially potent biomarkers to predict and track severity of disease. The importance of these novel cell populations is highlighted by their ability to robustly contribute to distinguishing severe and mild COVID-19 and acute COVID-19 from other viral infections (Fig. 4 and S7). Consequently, tracking these unique cells might provide important criteria for enrollment into clinical trials as well as provide a surrogate marker for tracking efficacy of new potential treatments. Finally, CPT1a has been associated with inflammasome activation, which has been observed in COVID-19 patients, and its role in fatty acid oxidation support modulators of either fatty acid oxidation or inflammasome signaling as a potentially novel therapeutic targets to mitigate disease. Together, our data demonstrate the utility of broad immuno-metabolic phenotyping to identify novel subsets of immune cells that have the potential to provide insight into disease pathogenesis and to define novel metabolic targets for treatment.

Funding

Funding for this work was provided by a Johns Hopkins University Provost Research Grant, The Bill and Melinda Gates Foundation (134582), NIH Centers of Excellence in Influenza Research and Surveillance (HHSN272201400007C) and NIH (P41EB028239 to JDP). The study was supported in part by a cooperative agreement between Johns Hopkins University (JHU) and the Division of Intramural Research, NIAID, NIH, as well as extramural support from the US Department of Health and Human Services Biomedical Advanced Research and Development Authority (BARDA; grant number IDSEP160031-01-00) and National Institute of Allergy and Infectious Diseases (R01AI120938, R01AI120938S1, R01AI128779, and U19A1088791). Contract HHSN272201400007C awarded to A.A.R.T the Johns Hopkins Center for Influenza Research and Surveillance (JHCEIRS) at the Johns Hopkins University. Confocal microscopy images acquired from Zeis LSM880 was supported by National Institute of General Medical Sciences (NIGMS) of the National Institutes of Health (S10OD023548 to JHUSOM Microscope Facility). The Core for Clinical Research Data Acquisition is supported in part by the Johns Hopkins Institute for Clinical and Translational Research (UL1TR001079).

Competing interests

J.D.P. is a scientific founder, a paid consultant and has equity in Dracen Pharmaceuticals.

Supplementary Materials

Materials and methods

Figures S1-S7

Tables S1-S2