Corticosteroid treatment in COVID-19 modulates host inflammatory responses and transcriptional signatures of immune dysregulation.

Severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) is the causative agent of coronavirus disease-2019 (COVID-19), a respiratory disease that varies in severity from mild to severe/fatal. Several risk factors for severe disease have been identified, notably age, male sex, and pre-existing conditions such as diabetes, obesity, and hypertension. Several advancements in clinical care have been achieved over the past year, including the use of corticosteroids (e.g., corticosteroids) and other immune-modulatory treatments that have now become standard of care for patients with acute severe COVID-19. While the understanding of the mechanisms that underlie increased disease severity with age has improved over the past few months, it remains incomplete. Furthermore, the molecular impact of corticosteroid treatment on host response to acute SARS-CoV-2 infection has not been investigated. In this study, a cross-sectional and longitudinal analysis of Ab, soluble immune mediators, and transcriptional responses in young (65 ≤ years) and aged (≥ 65 years) diabetic males with obesity hospitalized with acute severe COVID-19 was conducted. Additionally, the transcriptional profiles in samples obtained before and after corticosteroids became standard of care were compared. The analysis indicates that severe COVID-19 is characterized by robust Ab responses, heightened systemic inflammation, increased expression of genes related to inflammatory and pro-apoptotic processes, and reduced expression of those important for adaptive immunity regardless of age. In contrast, COVID-19 patients receiving steroids did not show high levels of systemic immune mediators and lacked transcriptional indicators of heightened inflammatory and apoptotic responses. Overall, these data suggest that inflammation and cell death are key drivers of severe COVID-19 pathogenesis in the absence of corticosteroid therapy. ©2021 Society for Leukocyte Biology.

Coronavirus Disease 2019

corticosteroids

coronavirus

days post symptom onset

Middle eastern respiratory syndrome

no corticosteroids treatment

receptor binding domain

severe acute respiratory syndrome

INTRODUCTION

Severe acute respiratory syndrome coronavirus 2 (SARS‐CoV‐2), the causative agent of coronavirus disease‐2019 (COVID‐19), is responsible for millions of deaths globally. Most infections result in asymptomatic or mild disease (∼80%). , , However, advanced age (>65 years) and pre‐existing comorbidities, such as obesity, diabetes, and pulmonary and cardio vascular disease, increase the probability of severe disease. , Severe COVID‐19 is characterized by pneumonia, lung fibrosis that often require ventilation, and a dysregulated inflammatory response that is followed by either recovery without complications, recovery with complications (“long COVID‐19″), or fatality. , Our understanding of COVID‐19 disease progression is continuously evolving as more studies investigate longitudinal changes in immune responses during severe COVID‐19. , ,

Infection with SARS‐CoV‐2 results in the sequential induction of the innate and adaptive immune responses. , Innate immunity following infection by the 4 endemic human coronaviruses (hCoVs) and closely related SARS‐CoV and Middle Eastern Respiratory Syndrome Coronavirus (MERS‐CoV) is induced by detection of viral RNA by pattern recognition receptors (PRRs) RIG‐1 and MDA‐5 and TLRs 3, 7, 8, and 9. , , , Recognition by these PRRs leads to production of type I IFN and transcription of IFN stimulated genes (ISGs) that are important for establishing an antiviral state. Importantly, SARS‐CoV‐2 is known to antagonize the type I IFN response, leading to dysregulated and altered immune responses. , Furthermore, a feature of severe COVID‐19 has been shown to be NLRP3 inflammasome activation leading to translation of IL‐1B, IL‐18, and pyroptosis, all of which enhance the pro‐inflammatory cytokine storm. , During the initiation of the adaptive immune response, antigen‐specific T and B cells are activated and recruited to the site of infection. , , The spike and nucleocapsid protein are the dominant targets of the antibody (Ab) response, while the main target of neutralizing antibodies is believed to be the receptor binding domain (RBD) of the spike protein. , In contrast, antigen‐specific T cells detect a variety of SARS‐CoV‐2 structural and nonstructural proteins. , , Multiple studies have performed correlative analyses that showed that severe disease is associated with heightened innate and dysregulated adaptive immune responses, however, the molecular signatures of these responses are unclear.

Several risk factors are associated with the development of severe COVID‐19. One of the most significant factors is age, with the overwhelming majority of hospitalizations and deaths occurring in individuals older than 65 years of age. , , , Additionally, obesity, diabetes, hypertension, and the male sex have also been associated with increased susceptibility to severe COVID‐19. , , , Earlier studies suggest that increased basal levels of systemic inflammation (“inflammaging”), reduced frequencies of naïve T and B cells, and dysregulated adaptive immune responses may facilitate disease progression. , , , For instance, a recent study found impaired cytotoxic CD8 T cell responses in aged, but not young patients, while others reported a positive correlation between age and inflammatory cytokines during acute COVID‐19. , However, the mechanisms underlying increased susceptibility to severe disease and dysregulated host responses to infection with age remain poorly defined.

In this study, we undertook a multiplatform approach to interrogate the host response in young (≤65 years) and aged (>65 years) males with severe COVID‐19 who also suffered from comorbidities known to increase susceptibility to severe disease (e.g., diabetes, obesity, cardiovascular diseases) during the acute phase of COVID‐19. Furthermore, we interrogated the impact of corticosteroid therapy, which has proven to be highly effective in suppressing the host inflammatory response. Specifically, we characterized humoral, immune mediators, and transcriptional responses using whole blood samples collected 1–3 weeks post symptoms onset. Additionally, RNA sequencing data were leveraged to infer changes in circulating immune cell populations. Collectively, the data presented here show dampened inflammatory and apoptotic responses in young and aged corticosteroid‐treated patients at the protein and transcriptional levels, suggesting that aberrant inflammation is a key driver of severe COVID‐19.

RESULTS

Severe COVID‐19 is associated with robust humoral responses and systemic inflammation that is modulated by corticosteroids

Fifty‐two blood samples were obtained from 41 male subjects categorized as “young” (≤65 years, n = 19) or “aged” (>65 years, n = 22) hospitalized with several COVID‐19. Diabetes, obesity, and cardiovascular disease were the most common comorbidities amongst our study participants (Supplementary Table S1). Patients whose samples were obtained between January and February 2021 had received corticosteroids, which had become standard of care by then (Graphical Abstract; Supplementary Table S1). To perform a longitudinal analysis, samples were grouped into “young” and “aged” groups and categorized into 2 different time (T) periods based on days post symptom onset (DPSO): T1 encompasses 6–11 DPSO, while T2 includes 15–21 DPSO. SARS‐CoV‐2‐specific antibodies (Abs) and soluble mediators were quantified longitudinally in plasma using ELISA and Luminex assay respectively, while the focus reduction neutralization test (FRNT) was used to measure neutralizing Ab titers against SARS‐CoV‐2 (Graphical Abstract; Supplementary Fig. S1A). Anti‐RBD IgG were readily detected at T1 in both age groups, with the modest increase at T2 in young patients only (Supplementary Fig. S1A). Neutralizing Ab titers showed no differences between age groups or timepoints (Supplementary Fig. S1B). Early in the pandemic, it has been reported that the bulk of the polyclonal neutralizing response targets the RBD region of SARS‐CoV‐2 spike protein and that RBD‐specific Ab titers correlate with neutralization. , , , , In line with these findings, we observed a positive correlation between anti‐RBD IgG endpoint titers and SARS‐CoV‐2 neutralization (r = 0.3902, P = 0.008) (Supplementary Fig. S1C). However, when we stratified patients by age, we observed a modest correlation between the RBD endpoint titers and live virus neutralization only in the young (r = 0.385, P = 0.070) compared to the aged cohort (r = 0.325, P = 0.140) suggesting that the polyclonal Ab response in the aged population differs from the adult population (Supplementary Fig. S1C).

We next measured levels of circulating immune mediators since severe COVID‐19 is associated with a cytokine storm (Supplementary Fig. S1D and E). Additionally, patients were divided into two groups based on corticosteroid treatment (CORT) treatment or not (NC), and immune mediator levels were compared to those measured in age‐matched male healthy donors (HDs). Samples obtained before corticosteroid became standard of care exhibited significantly elevated levels of numerous cytokines, chemokines, and growth factors compared to HDs (Supplementary Fig. S1D). Levels of both pro‐inflammatory (e.g., IL‐1α/b, IL‐6, IL‐17a, IL‐33) and anti‐inflammatory (e.g., IL‐13) cytokines were significantly increased in both age groups (Supplementary Fig. S1D). Additionally, cytokines important for leukocyte proliferation and survival (e.g., IL‐2, IL‐3, IL‐7, IL‐12p70) were higher in both young and aged patients (Supplementary Fig. S1D). Markers of T cell‐mediated immunity including CD40L, granzyme B, and PDL1 were elevated in young and aged patients, with CD40L produced to a significantly greater level in young compared to aged patients. Although there was no significant increase in levels of interferon‐α (IFNα) relative to HD, we observed enhanced levels of IFN‐β and IFN‐γ in both age groups (Supplementary Fig. S1D). Chemokines primarily involved in the recruitment of monocytes (e.g., fractalkine, G‐CSF, MIP‐1α) were elevated regardless of age, while MIP1b was exclusively elevated in aged patients. Mitogenic growth factors, such as PDGF‐AA and VEGF, were notably elevated in all patients (Supplementary Fig. S1D). In contrast, fewer immune mediators were increased in CORT‐treated young and aged patients (Supplementary Fig. S1E). We detected significantly higher concentrations of IL10, TNFα, IP10, RANTES, and EGF in both young and aged CORT‐treated patients compared to HD (Supplementary Fig. S1E). Two cytokines involved in anti‐inflammatory and cellular cytotoxic responses, IL1RA and IL15 respectively, were significantly elevated in young CORT‐treated patients only. On the other hand, levels of the chemokine, MIP1β were significantly higher in aged CORT‐treated patients only (Supplementary Fig. S1E). Unlike NC patients, no significant changes in levels of IFN‐γ and IFN‐β were noted in CORT‐treated young and aged patients.

Steroid treatment modulates host transcriptional responses during acute COVID‐19

We next investigated the host transcriptional response to severe COVID‐19 in both age groups. Transcriptional profiles of the patients were compared to that of sex‐ and age‐matched healthy donors (HD, NCBI BioProject PRJNA511612). Principal component analysis (PCA) showed a distinct separation between HD and COVID‐19 patients, as well as between CORT and NC patients regardless of age, with no clear distinction between T1 and T2 (Supplementary Fig. S2A).

We first analyzed the transcriptional response in NC patients compared to age‐matched HD (Fig. 1 ; Supplementary Fig. ). Acute SARS‐CoV‐2 infection results in substantial gene expression changes with ∼2000 differentiatlly expressed genes (DEGs) (1344 upregulated; 933 downregulated in young ND; and 1397 upregulated; 569 downregulated in aged ND), with a large overlap between the 2 age groups (Fig. 1A). Functional enrichment of DEGs expressed by both groups showed over‐representation of Gene Ontology (GO) terms characteristic of antiviral host defense, including “viral gene expression” (e.g., RPL24, RPS7), “myeloid leukocyte activation” (e.g., CCR2, CD177, MMMP8/9, CXCL1/8, S100A8/11/12), “lymphocyte activation” (e.g., BCL6, CD28, CD3D, IGHA1/A2/G1/G2/G3/4, LYN), and “inflammatory response” (e.g., IL17RA, IL1B, RELB, RIPK2, TGFB1) (Fig. 1B; Supplementary Fig. S2A and B). Only DEGs from young patients enriched to GO term “antigen processing and presentation,” which included a number of genes important for MHC class I/II presentation (e.g., CD1D/74, HLA‐G/DR/DOA/DMB/DQA1//DQB1) and peptide processing (e.g., TAP2, TAPBP/L, ERAP1/2, CTSD, PSMA5) (Fig. 1B and C). In contrast, only DEGs detected in the aged patients enriched to GO terms associated with adaptive immunity (e.g., “B cell mediated immunity,” “CD4 αβ T cell activation”). These DEGs encoded proteins with roles in T cell differentiation (e.g., IL4R, JAK3, RORA, ZBTB7B), activation (e.g., ARS, CD86), and regulation (e.g., FOXP1, TCIRG1), and were primarily upregulated relative to HD (Fig. 1D).

FIGURE 1

Transcriptional profiling of severe COVID‐19 reflects innate and adaptive immune dysregulation. (A) Venn diagram of differentially expressed genes (DEGs) computed for young or aged patients not treated with corticosteroids (ND) relative to age‐ and sex‐matched HD. (B) Bubbleplot representing the functional enrichment of DEGs detected in young and aged patients identified in panel (A). The size of each circle represents the number of DEGs belonging to the indicated gene ontology (GO) term while color represents the FDR‐corrected P‐value (q‐value). (C) Violin plots depicting normalized transcript counts of DEGs belonging to GO term “antigen processing and presentation” (young patients); grey color represents healthy donors (HD) while the blue color represents the non‐CORT treated group (NC) (D) Violin plots depicting normalized transcript counts of DEGs belonging to GO terms “CD4 alpha‐beta T cell activation” and “B cell‐mediated immunity” (aged patients); grey color represents healthy donors (HD) while the blue color represents the non‐CORT treated group (NC)

Next, we carried out the same analysis for the larger cohort of CORT‐treated patients to examine the effects of CORT on transcriptional changes (Supplementary Fig. 2; Supplementary Fig. ). Comparisons of each CORT‐treated age group to age‐matched HD identified nearly 1500 DEGs in each group (Fig. 2A). The majority of DEGs were upregulated in both young (84%) and aged (89%) patients with a significant overlap but were not differentially expressed between the 2 age groups (Fig. 2A). DEGs expressed in young and aged patients enriched to GO terms primarily associated with antiviral innate immunity (e.g., “response to virus”; DDX58, EIF2AK2A, IFI16, IRF7, MX1, OAS2/3, STAT2); inflammation (“myeloid leukocyte activation”; C5AR1, CCR1, CD177, MMP8/9/25, MPO, S100A9/11); and cell death (“positive regulation of programmed cell death”; ADAM8, BCL2L1, GZMA) (Fig. 2B,D; Supplementary Fig. S3A). Additionally, DEGs detected in both age groups enriched to GO terms related to adaptive immunity, such as “B cell mediated immunity” (e.g., IGHA, IGHM, IGHG1/2/3/4, IGHV3‐21, IGKC, IGLC3) and “lymphocyte activation” (e.g., ARG1, BCL6, CCR7, CD3D, HLA‐B, LCK, JAK3; Fig. 2C; Supplementary Fig. S3B). DEGs detected only in young patients enriched to GO terms related to inflammation such as “interleukin‐6 production” and “wound healing” (e.g., CCR2, CX3CR1, IL17RA, TGFB1, HBB, ITGA5, PLAUR, THBS1; Supplementary Fig. S3C). While DEGs detected only in aged patients also enriched to inflammation‐related GO terms (e.g. “I‐kappaB kinase/NF‐kappaB signaling”; BCL3, NFKBIA, NFKBID, RELB) as well as adaptive immunity (“antimicrobial humoral response”, “CD4 αβ T cell differentiation”; CD83, SATSH3, ZBTB7B; Fig. 2C ; Supplementary S3D).

FIGURE 2

Young and aged corticosteroid‐treated COVID‐19 patients exhibit overlapping and distinct transcriptomes. (A) Venn diagram of differentially expressed genes (DEGs) computed for CORT‐treated patients relative to age‐ and sex‐matched HD in panel (A). (B) Bubbleplot representing the functional enrichment of DEGs detected in young and aged patients identified in panel (B). The size of each circle represents the number of DEGs belonging to the indicated gene ontology (GO) term while color represents the FDR‐corrected P‐value (q‐value). (C) Heatmap representing DEGs enriching to all DEGs enriching to “response to virus” in young and aged patients. Each column represents the average rpkm for the given group. Range of colors is based on row‐scaled and centered rpkm values of the represented DEGs. Red represents high expression while blue represents low expression

To further investigate the impact of CORT treatment on COVID‐19 progression at the molecular level, we compared the transcriptional profiles between CORT‐treated and NC patients (Fig. 3). Significant overlap was detected between the transcriptional response to SARS‐CoV‐2 infection generated by untreated and CORT‐treated patients within each age group (Fig. 3A and B). DEGs unique to or shared between the treatment groups enriched to similar GO terms indicative of robust innate (e.g. “myeloid leukocyte activation/differentiation”, “response to bacterium”) and adaptive (e.g., “T cell differentiation”) immunity (Fig. 3C and D). For instance, in young patients (regardless of treatment) expression of genes important for T cell immunity such as CD28, and HLA‐DMB were downregulated while that of B cell genes (e.g., IGHA1/2, IGHG1/2/3/4) and pro‐apoptotic genes (e.g., BCL2A1, CD177) were upregulated ( ). In untreated young patients, other genes important for antiviral defense were downregulated (e.g., CD27, CD3G, DDX58), while others were upregulated in CORT‐treated patients (e.g., HLA‐B/E/F, IL4R, RSAD2; Supplementary Fig. S4A). A subset of apoptotic genes, notably CASP8 and FADD, were highly induced only in untreated young patients. Aged NC and CORT‐treated patients also shared a core of highly expressed lymphocyte activation genes such as FOS and IGHG1/2/3/4 (Supplementary Fig. S4A). A subset of T‐ and B‐cell genes (e.g., HLA‐A/B/E/F, IGHV3‐21) as well as inflammatory genes such as IL18RA and S100A9/P associated with neutrophil‐mediated immunity were uniquely upregulated in CORT‐treated aged patients (Supplementary Fig. S4B). Genes encoding pro‐inflammatory factors and ISGs, including BCL2L11, BIRC2, IFIH1, and IFNGR2, were exclusively expressed to high levels in only aged NC patients.

FIGURE 3

Corticosteroids treatment modulates transcriptional responses in acute COVID‐19. Venn diagrams representing DEGs between corticosteroid (CORT)‐treated and untreated (NC) for (A) young and (B) aged patients. Functional enrichment of CORT‐unique, NC‐unique and shared DEGs in panels (A) and (B) for (C) young and (D) aged cohorts. Color intensity represents the statistical significance (shown as the negative log of the FDR‐adjust P‐value, ‐log(q‐value)) with a range of colors based on the GO terms with the lowest and highest –log(q‐value) values for the entire set of genes per cluster. Numbers of genes enriching to each GO term per cluster are represented in each box; blank boxes indicate no statistical significance

Longitudinal regression analyses identify genes involved in antiviral and immune processes

In order to identify clusters of genes that are similarly regulated overtime from T1 to T2 in both young and aged CORT‐treated patients, we employed a two‐way forward stepwise regression analysis (Fig. 4). This analysis was only carried out using samples from CORT‐treated patients since we had a large number of samples in this category. Genes considered significant in a minimum of 16 comparisons were retained and clustered by temporal expression patterns, which resulted in two clusters (Fig. 4A). The 768 genes in cluster 1 were downregulated at T1 and their expression remained reduced from T1 to T2, while genes in cluster 2 (n = 73) displayed the opposite pattern (Fig. 4A). Genes in cluster 1 enriched to several terms illustrating lymphocyte‐mediated immunity (e.g., “lymphocyte differentiation,” “TCR signaling”) and response to signaling/stress (e.g., “response to peptide,” “macroautophagy”; Fig. 4B). Likewise, we detected several genes encoding subunits of the T cell receptor (e.g., TRAJ12, TRAC) as well as T cell coreceptors like CD28 and CD247 (Fig. 4C). Several apoptosis mediators (e.g., TNFRSF10B) and genes regulating B cell immunity (e.g., MEF2C) also belonged to cluster 2. The few genes in cluster 2 included growth factor TGFB1, protease FURIN, and NFKB transcription factor component RELA, all of which were expressed to similar levels by young and aged patients at T1 and T2 (Fig. 4D).

FIGURE 4

Regression analysis identifies longitudinal downregulation of key antiviral and T cell genes in both young and aged CORT‐treated patients. (A) Gene expression (in rpkm) of 2 gene clusters identified by 2‐ways forward regression analysis. (B) Gene ontology (GO) network depicting functional enrichment of the 500 genes contributing the most to variation seen in PCA using Metascape. Clustered nodes of identical color correspond to one GO term. Node size represents the number of DEGs associated with the GO term. Gray lines represent shared interactions between GO terms, with density and number indicating the strengths of connections between closely related GO terms. (C) Heatmap depicting gene expression of cluster 1 genes. Exemplar genes are labeled. Each column represents the average rpkm for the given group and timepoint. Range of colors is based on row‐scaled and centered rpkm values of the all the DEGs in the heatmap. Red represents increased while blue represents decreased expression relative to the median. (D) Beanplot representing expression of exemplar genes identified in cluster 2 (panel A) in young and aged patients at T1 and T2

Next, we refined our longitudinal analysis by performing a regression against DPSO to identify longitudinally regulated genes (defined as protein‐coding genes with FDR ≤ 0.05; Fig. 5). We first carried out this regression independently for each age group (Fig. 5A and B). This analysis uncovered 108 genes in the CORT‐treated aged group only. These genes were mostly downregulated overtime and enriched to GO terms related to classical complement (e.g., IGHV1‐18, IGHV2‐70), TLR signaling pathways (e.g., BIRC2, IRF4), as well as epigenetic modification and cell division (e.g., DNAJC6, ESCO1, H2AC12) (Fig. 5A and B). We then carried out these analyses using age as a co‐variate. We identified 102 genes, 3 of which, DDX6 (viral RNA sensing), RAB2A (intracellular trafficking), and UBE2S (modulation of viral replication), showed opposite expression trends in the 2 age groups (Fig. 5C). VAC14, which encodes a scaffold protein for intracellular vesicle transport, was downregulated in both young and aged patients over time, albeit to a greater magnitude in aged patients (Fig. 5C).

FIGURE 5

Genes related to humoral immunity and DNA metabolism are longitudinally downregulated in aged CORT‐treated patients. (A) Functional enrichment of LRGs identified for aged CORT‐treated patients only. Color of each bar color represents the FDR‐corrected P‐value (q‐value). (B) Expression of LRGs from panel (B). Each column represents the average (or singular value) rpkm for the given timepoint. Range of colors is based on row‐scaled and centered rpkm values of the represented DEGs. Red represents upregulation; blue represents down‐regulation. (C) Expression of longitudinal‐regulated genes (LRGs) (DD6, RAB2A, UBE2S, VAC14) longitudinally regulated in aged and young corticosteroids (CORT)‐treated patients

Severe COVID‐19 is associated with profound immune cell dysregulation

Since these blood samples were obtained from a biobank premixed with an RNA stabilizer, immune phenotyping of circulating immune cells could not be performed in this study. Therefore, we applied deconvolutional analysis to our bulk RNA sequencing (RNA‐seq) data to infer changes in cell frequencies (Supplementary Fig. S5). To validate this approach, we obtained bulk RNA‐seq and single‐cell (sc)RNA‐seq from a previously published study. We used scRNA‐Seq cell type‐specific marker genes to deconvolute the bulk RNAseq data from the same cohort. The deconvoluted cell type proportions matched the proportions in the single‐cell data set remarkably well, giving credence to the validity of the deconvolution procedure (Supplementary Fig. S5A). , , , We first compared cell frequencies among HD and CORT‐treated patients. Although no significant differences were observed, there was a trend towards a decrease in CD4 memory cells and an increase in CD14 monocytes in severe COVID‐19 patients irrespective of their age (Supplementary Fig. S5B and Table S3). No significant differences in the proportions of major immune cell types were detected between young and aged patients other than a modest decrease in CD8 naïve lymphocytes and CD16+ monocytes in aged patients (Supplementary Fig. S5C and Table S4).

To further delineate differences in cell subsets between HD and COVID‐19 patients at a higher resolution, we next performed digital cell quantification (DCQ) using in silico flow cytometry. We used the bulk RNA‐Seq data to infer cell proportions using ImmQuant with the IRIS database (Supplementary Fig. S5D and E ). , Our analysis revealed a dysregulation of both innate and adaptive immune cell subsets in young and aged patients compared to HD. Specifically, significant increases in bone marrow‐derived plasma cells were predicted for all patients regardless of treatment (Supplementary Fig. S5D and E). We also observed a significant decrease in IL15‐stimulated NK cells, activated day 7 monocytes, and CD4 Th1 cells, whereas there were significant increases in CD4 Th2 and plasma cells, young and aged NC patients only. Increases in monocytes and neutrophils were anticipated at later timepoints in CORT‐treated patients. A unique decrease in IL2‐stimulated NK cells was predicted for young CORT‐treated patients at T1‐T2 (Supplementary Fig. S5D and E). No significant differences were observed for cell frequencies between T1 and T2.

DISCUSSION

Our current understanding of COVID‐19 is rapidly evolving. However, immune mechanisms governing severe disease remain elusive but critical for treating short‐ and long‐term consequences of COVID‐19 in vulnerable populations. In this study, we performed a longitudinal analysis of the functional and transcriptional host responses to severe COVID‐19 in the peripheral blood of young (≤65 years) and aged (≥65 years) male patients who suffered from co‐morbidities known to exacerbate COVID‐19 severity. Furthermore, to our knowledge, this is the first study to investigate the effects of corticosteroid therapy (i.e., corticosteroids; CORT) on host responses during severe COVID‐19.

Virus‐specific Abs are critical determinants of disease progression and overall immunity. Multiple studies reported a strong induction of IgG against RBD of SARS‐CoV‐2 spike protein within 2 weeks post symptom onset, where anti‐RBD and anti‐spike IgG titer magnitudes are directly correlated with disease severity. , , , , Indeed, we observed an induction of RBD‐specific IgG Abs within the first 6–11 days post symptom onset (DPSO) in both young and aged NC and CORT‐treated patients that correlated with neutralizing Abs (nAbs) titers, primarily in young patients. These data support studies suggesting that many neutralizing Abs are directed toward the RBD of the spike to prevent viral entry. , , , , However, in the aged cohort there is no clear correlation between the level of RBD binding Abs and neutralization. This could be occurring due to an increase in non‐neutralizing RBD‐specific Abs or an increase in neutralizing Abs outside of RBD within the aged cohort. Interestingly, the development of antigen‐specific IgG and neutralizing Abs occurred simultaneously with predicted increases in the frequencies of plasma B cells in untreated patients while upregulation of B cell‐related genes was most notable in CORT‐treated patients. The robust transcriptional inflammatory response accompanied by high levels of systemic immune mediators, and a significant downregulation of immunoglobulin genes by NC patients, suggest a potential for poor B cell responses and decreased Ab durability at later time points. , Poor Ab durability in severe disease may be due to dysregulated germinal center responses driven by a milieu of inflammatory cytokines and extra‐follicular B cell activation that generates short‐lived plasma cells. , We observed significantly elevated levels of TNF‐α, a cytokine known to drive germinal center dysregulation, in all patients regardless of age or treatment and preliminary studies suggested beneficial effects of anti‐TNF‐α therapy. , , Protection offered by peripheral Abs is still unclear, but a correlation between plasma and mucosal Abs suggests that protective Ab responses also be detected at the site of infection.

Consistent with previous studies, we detected a robust cytokine storm in severe COVID‐19 patients that was significantly attenuated by corticosteroid treatment. , , , , , , , , , , , , , In untreated patients, elevated levels of cytokines previously associated with severe COVID‐19 were noted, including IL‐1b, IL‐6, IL‐10, IL‐12, IFN‐γ, and TNF‐α as previously reported. , , , , , , , , Moreover, levels of TNF‐α and MIP3β were higher in aged patients relative to young patients, suggesting that some cytokines may serve as early, age‐specific prognostic markers of severe disease. Indeed G‐CSF and GM‐CSF, which were significantly increased in NC patients, were previously identified as hallmarks of severe COVID‐19 in other studies. , , , , In contrast, a limited number of pro‐ and anti‐inflammatory cytokines, chemokines, and growth factors were detected in CORT‐treated patients. Corticosteroids became standard of care at the end of June 2020 given the significant reduction in 28‐day mortality and the need for mechanical ventilation they provide. In CORT‐treated patients, levels of several inflamamtory cytokines and chemokines were elevated (TNF‐α, eotaxin, IP‐10, RANTES, and EGF). The anti‐inflammatory IL‐1RA, NK cell/memory T cell survival factor IL‐15, and monocyte attractant MCP‐1 were elevated in young patients only while MIP‐1β, which is involved in the recruitment of NK and T cells, was elevated in CORT‐treated aged patients. These observations indicate that despite comparable disease severity scores, age impacts the nature and magnitude of the cytokine response to infection. It remains unclear whether peripheral cytokines and chemokines originate from lung‐resident versus peripheral blood cells.

Our bivariate and regression transcriptional analyses of whole blood samples from CORT‐treated patients revealed robust transcriptional changes indicative of dysregulated inflammatory responses despite reduced levels of systemic inflammation. Moreover, in contrast to the large transcriptional differences observed in blood samples obtained earlier in the pandemic before CORT treatment became standard of care, the gene expression profiles of young and aged samples were more uniform post CORT treatment. Genes encoding proteins that play a role in coagulation and wounding were also upregulated in CORT‐treated patients, suggesting that long‐term complications including thrombotic events may not be reduced with corticoid treatment. Genes encoding members of the NFKB signaling pathway were largely induced in the aged cohort (e.g., NFKBIA), thrombotic pathways in the young (e.g., THBS1), and complement pathway in both groups (e.g., C5AR1). , , , , These gene expression changes are predicted to parallel increased frequency of neutrophils and monocytes, which are a major sources of inflammatory mediators in COVID‐19. , Indeed, genes associated with neutrophil defense and extracellular traps (NETs) (e.g., MPO, S100A12) were highly upregulated in both CORT‐treated age groups. Severe COVID‐19, as well as other severe viral respiratory infections, have been associated with NET formation, which may lead to enhanced inflammation, pulmonary dysfunction and changes in vasculature. , , , , , Interestingly, we did not detect significant increases in systemic levels of IL‐8, a neutrophil‐specific chemoattractant and inducer of NET formation.

A deficient antiviral IFN response has been proposed as the mechanism behind severe COVID‐19. , , , Levels of type I (IFN‐β) and II (IFN‐γ) IFN, but not IFN‐α, were highly elevated in NC patients and accompanied by increases in cytokines that drive IFN production (e.g., IL‐2, IL‐12, and IL‐13) and activation of IFN‐secreting immune cells (e.g., IL‐3, IL‐7, IL‐33). In contrast, levels of IFN‐β and IFN‐γ, were not significantly elevated compared to healthy donors in CORT‐treated patients. However, the expression of IFN‐stimulated genes (ISGs) was comparable between CORT‐treated and untreated patients with only a few ISGs downregulated in untreated patients, such as STAT2 and STING1 recently implicated in RNA virus detection. These data suggest that corticosteroid treatment can modulate type I IFN production, but not enough to disrupt ISG expression, which is likely initiated early in the infection (before CORT treatment). IFN production in NC patients could be mediated by peripheral plasmacytoid dendritic cells (pDCs; major producers of types I IFNs) and T cells (major producers of type II IFNs) as has been reported by other studies. , , , , , , , ,

On the other hand, we detected robust changes in the expression of T cell activation genes in CORT‐treated patients, suggesting dysregulation in cellular immunity and in line with our predicted T cell lymphopenia as reported by other studies. , , , , , , , , We also detected lower expression of genes critical for T cell‐mediated immunity in NC patients including a number of MHC class I molecules (e.g., HLA‐A) and T cell receptor signaling components (e.g., ZAP70). Additionally, we detected significantly elevated levels of soluble PDL1 in the plasma of NC patients, suggesting a dysregulation of adaptive immunity and potential exhaustion in these patients. The role of the antiviral response in the periphery is still unclear given that SARS‐CoV‐2 replication is restricted to the respiratory tract and that actively replicating virus has not been detected in the blood.

There are several limitations to our study. First, our study strictly examined severe COVID‐19 in male individuals ages 46 to 76 with obesity and diabetes. Given the sex‐specific differences in the immune responses to SARS‐CoV‐2 infection, we would therefore, caution against generalization of the data reported in this study. Second, given our focus on only severe COVID‐19 patients, we are unable to identify statistically significant biomarkers that differentiate mild, moderate, and severe COVID‐19 cases from each other. Furthermore, we lack a comparison to COVID‐19 patients without co‐morbidities to uncover mechanisms driven by pre‐existing medical conditions, although such a cohort is unlikely given that obesity, diabetes and hypertension are significant risk factors for severe COVID‐19. , Third, our findings specifically examine the peripheral blood response in severe COVID‐19 rather than the site of infection: the respiratory tract. Also, our study does not profile post‐acute phase of COVID‐19 to understand short‐ and long‐term outcomes (e.g., “long COVID19”). Future studies should aim to compare longitudinal respiratory and peripheral samples collected from larger cohorts of individuals with varying disease severity to identify and characterize tissue‐specific transcriptional and functional changes. This will enable the identification of biomarkers of different disease severities. Our study also provide data to inform future clinical trials of COVID‐19 therapeutics and advocates for the use of corticosteroids, such as dexamethasone. We observed that levels of IL6, a prognostic marker for the development of severe COVID‐19 and fatality, were unchanged in CORT‐treated, but not NC patients. The long‐term impact of CORT treatment on long‐term immunity and vaccine responsiveness needs further investigation.

In summary, we performed a multi‐platformed study incorporating both functional and transcriptional data to provide crucial insight into the molecular pathogenesis of severe COVID‐19 in the peripheral blood of patients receiving and not receiving corticosteroids treatment. Our data support the hypothesis that severe COVID‐19 is associated with pronounced dysregulation of all major immune cell subsets in blood, which confirms immunological misfiring in severe COVID‐19. However, profound changes are mediated by corticosteroid therapy. The synergy of dysregulated innate and adaptive immunity fueled by an aberrant cytokine response may be key to driving systemic pathology in severe COVID‐19.

MATERIALS AND METHODS

Study design

Remnant blood samples (n = 36) were collected from hospitalized male patients with severe COVID‐19 from 2020 to 2021 through the COVID biobank of the university of California Irvine. COVID‐19 diagnosis was confirmed positive with a nasopharyngeal positive qRT‐PCR test. Samples were classified into 2 time points based on days post symptom onset (DPSO): T1 (6‐11 DPSO) and T2 (15‐21 DPSO) (Fig. 1A). Samples were further classified into young (≤65 years) and aged (≥65 years) cohorts based on age. Complete metadata for COVD19 patients and healthy donors can be found in Supplementary . This study was approved by University of California Irvine Institutional Review Board (HS# 2012–8716).

ELISA endpoint titer

SARS‐CoV‐2 specific IgG titers against the receptor‐binding domain (RBD) of the spike protein was determined by ELISA. For this, MaxiSorp (ThermoFisher) plates were coated with 50 μL of a 1 μg/mL mixture of protein in carbonate buffer (0.1 M Na2CO3 0.1 M NaHCO3 pH 9.3) overnight at 4°C. Blocking buffer (PBS + 5%BSA + 0.5% Tween) was then added to pates for 1 h at 37°C followed by washing 3× with wash buffer. Serially diluted polyclonal sera were then added at a volume of 50uL/well. Polyclonal sera were incubated for 1 h at room temperature in the ELISA plate, washed 3x with wash buffer, followed by addition of goat‐anti‐human IgG HRP (Sigma) conjugated secondary (1:5000) for 1 h at room temperature. The plate was washed 3× with wash buffer and the ELISA was developed with 100 μL of TMB enhanced substrate (Neogen Diagnostics) and placed in a dark space for 15 min. The reaction was quenched with 1 N HCl and the plate was read for an optical density (OD) of 450 nm on a BioTek Epoch plate reader. IgG endpoint titer were calculated using GrapPad Prism 8.

Focus reduction neutralization test

Human sera were diluted 4‐fold in a serial manner and mixed with ∼100 focus‐forming units (FFU) of virus followed by incubation at 37°C for 1 h to allow for the formation of immune complexes. The sera and virus mixture were then added onto confluent Vero‐WHO monolayers in 96‐well plates for 1 h at 37°C to allow virus adsorption. Cells were then overlaid with 2% methylcellulose mixed with DMEM containing 5% FBS and incubated for 24 h at 37°C. After 24 h, the media was removed, and the monolayers were fixed with 5% paraformaldehyde in PBS for 15 min at room temperature, rinsed with PBS, and permeabilized in Perm Wash (PBS, 0.05% Triton‐X). Infected cell foci were stained by incubating with polyclonal anti‐SARS Guinea Pig sera (BEIresources) for 1 h at room temperature and then washed 3 times with Perm Wash. Detection of foci was completed by incubation with a 1:5000 dilution of horseradish peroxidase‐conjugated goat anti‐guinea pig IgG (Sigma) for 1 h. After 3 washes with Perm Wash, staining was visualized by addition of TrueBlue detection reagent (KPL). Infected foci were enumerated by CTL Elispot and FRNT curves were generated by log‐transformation of the x‐axis followed by non‐linear curve fit regression analysis using Graphpad Prism 8.

Luminex assay

Soluble immune mediators in plasma were quantified with Human XL Cytokine Luminex Performance Panel Kit (Thermo Scientific) per manufacturer's instructions. This panel quantified levels of the following soluble mediators: IL1ra, IL1b, IL1a, IL2, IL3, IL33, IL4, IL5, IL6, IL7, IL8, IL10, IL12p70, IL13, IL15, IL17A, IL17E, RANTES (i.e., CCL5), TNF‐α, IFN‐α, IFN‐β, IFN‐γ, eotaxin, fractalkine, G‐CSF, GM‐CSF, IFN‐γ‐induced protein 10 (IP10), Mϕ chemoattractant protein‐1 (MCP1), Mϕ inflammatory protein 1α (MIP1α), Mϕ inflammatory protein 1β (MIP1β), Mϕ inflammatory protein 3α (MIP3α), growth‐regulated oncogene α (GROα), growth‐regulated oncogene β (GROβ), platelet‐derived growth factor A (PDGF‐A), platelet‐derived growth factor AB/BB (PDGF‐AB/BB), transforming growth factor α (TGF‐α), tumor necrosis family‐related apoptosis‐inducing ligand (TRAIL), vascular endothelial growth factor (VEGF), endothelial growth factor (EGF), basic fibroblast growth factor (FGF‐basic), Fms‐related tyrosine kinase 3 ligand (Flt‐3 L), CD40L, and granzyme B. Analytes with undetectable levels in the majority of patients were omitted from analysis. COVID‐19 patients were compared to age‐matched healthy male donors.

Library construction and sequencing

RNA was extracted from whole blood using the RiboPure™‐Blood Kit (Life Technologies) per the manufacturer's instructions. Quality and quantity of RNA from whole blood was determined using an Agilent 2100 Bioanalyzer. cDNA libraries were constructed using the NEB Next Ultra II Direction RNA Library Prep Kit (Thermo Fischer). RNA was treated with RNase H and DNase I after depletion of ribosomal RNA (rRNA). Adapters were ligated to cDNA products and the subsequent ∼300 base pair (bp) amplicon was PCR‐amplified and selected by size exclusion. cDNA libraries were assessed for quality and quantity on an Agilent 2100 Bioanalyzer prior to 150 bp single‐end sequencing using the Illumina NovaSeq platform.

Bioinformatic analysis

Preliminary data analysis was performed with RNA‐Seq workflow module of systemPipeR. RNA‐Seq reads were demultiplexed, quality‐filtered and trimmed using Trim Galore (average Phred score cut‐off of 30, minimum length of 50 bp). FastQC was used to generate quality reports. Hisat2 was used to align reads to the reference genome Homo sapiens (Homo_sapiens.GRCh38.dna.primary_assembly.fa) and the Homo_sapiens.GRCh38.85.gtf file was used for annotation. Raw expression values (gene‐level read counts) were generated using the summarizeOverlaps function and normalized (read per kilobase of transcript per million mapped reads, rpkm) using the edgeR (v.3.30.3) and systemPipeR (v1.22.0) packages. Differentially expressed genes (DEGs) were defined as protein‐coding genes with an average rpkm ≥ 5, FDR‐corrected P‐value ≤ 0.05 and log2(fold‐change) ≥ 1 in either direction. To identify significant patterns of longitudinal gene expression, we applied a two‐ways forward regression analysis using MaSigPro. Only clusters containing genes that were considered statistically significant in at least 16 comparisons were retained. as the total number of normalized reads mapping to all viral genes.

Longitudinal analysis for CORT‐treated young and age patients over T1 and T2 was performed with short time‐series expression miner (STEM). STEM specifically integrates a clustering algorithm based on the sole assumption that experiments can be naturally and sequentially ordered. This software further differentiates between random and chronological patterns. Identification of longitudinally regulated genes (LRGs) with age as a factor or by age cohort in CORT‐treated patients was performed with DESeq2.

Functional enrichment of genes was performed with Metascape to identify significant Gene Ontology (GO) biological processes. The Cytoscape network data integration and visualization tool was used to generate GO (biological process) term network. Bar graphs, beanplots, bubbleplots, heatmaps, Venn diagrams and volcano plots were created in R. Graphs were generated with GraphPad Prism software (version 8).

Digital cell quantification (DCQ), (i.e., in silico flow cytometry) using gene expression data was carried out with ImmQuant using the IRIS database. Transcriptional data for age‐ and sex‐matched healthy donors were retrieved from NBCI BioProject PRJNA511612.

Deconvolution of RNA‐Seq data

In order to estimate the proportion of major cell subsets in blood samples, cell type‐specific marker genes were identified using a reference scRNA‐Seq dataset and used to deconvolute bulk RNAseq expression data. The scRNAseq data and cluster annotations were obtained from a recent study; bulk RNAseq data were generated in the present study. Both the bulk and single cell count data were normalized using the vst function from the sctransform R package. Cell type‐specific markers for each scRNAseq cluster that were used for deconvolutional analysis were found using the NS‐Forest machine learning method. Marker gene by sample bulk and marker gene by cell single‐cell matrices were used for deconvolution using the damped weighted least squares method.

Statistical analysis

Statistical significance between groups in panels 1A, 1B, 1D and 1E were computed with Mann‐Whitney U‐test panels S5A and S5B with one‐way ANOVA with Tukey's HSD. Spearman correlation was used to calculate correlations in panel 1C. *P ≤ 0.05, ** P ≤ 0.01, *** P ≤ 0.001, **** P ≤ 0.0001.

AUTHORSHIP

I.M., A.N.P., A.K.P., and J.D.B. conceptualized the work; T.L.S., J.D.B., A.N.P., K.J.M., and M.Z.Z. carried out the experiments; T.L.S., J.D.B., A.N.P., M.Z.Z., A.K.P., Y.Z., N.T., R.H.S., and I.M. analyzed the data; T.L.S., A.K.P., J.D.B., A.N.P., M.Z.Z., and I.M. wrote and edited the manuscript. A.D., S.L.G., and D.T. were responsible for the collection of patient samples. S.L.G. is the PI for the ACTTI trial. All authors reviewed and approved the final version of the manuscript.

DISCLOSURE

The authors declare no competing interest.

AVAILABILITY

The accession numbers for the RNA‐Sequencing data reported in this paper is available in the NCBI BioProject PRJNA751085 .

Figure S1. Corticosteroids attenuates the cytokine storm in severe COVID‐19.

Click here for additional data file.

Figure S2. Transcriptional analysis of acute COVID‐19 in young and aged patients.

Click here for additional data file.

Figure S3. CORT‐treated young COVID‐19 patients feature unique upregulation of genes modulating coagulation, thrombosis and inflammatory responses.

Click here for additional data file.

Figure S4. ND COVID‐19 patients are associated with unique upregulation of inflammatory and pro‐apoptotic genes.

Click here for additional data file.

Figure S5. Digital cell quantification predicts severe lymphopenia.

Click here for additional data file.

Table S1. Participant demographics and clinical data.

Table S2. Concentrations (pg/mL) of soluble immune mediators in young and aged COVID‐19 patients and healthy donors. NA, not detected or quantifiable.

Table S3. Statistical analysis of deconvolution results between HD and CORT‐treated COVID19 patients.

Table S4. Statistical analysis of deconvolution results between young and aged CORT‐treated COVID19 patients.

Click here for additional data file.