Abstract
We performed comparative lower respiratory tract transcriptional profiling of 52 critically ill patients with ARDS from COVID-19 or other etiologies, or without ARDS. We found no evidence of cytokine storm but instead observed complex host response dysregulation driven by genes with non-canonical roles in inflammation and immunity that were predicted to be modulated by dexamethasone. Compared to other viral ARDS, COVID-19 was characterized by impaired interferon-stimulated gene expression.
Brief Communication
In its most severe form, coronavirus disease 2019 (COVID-19) can precipitate the acute respiratory distress syndrome (ARDS), which is characterized by low arterial oxygen concentrations, alveolar injury and a dysregulated inflammatory response in the lungs1. Early reports hypothesized that COVID-19 ARDS was driven by a “cytokine storm” based on the detection of higher circulating inflammatory cytokine levels in critically ill COVID-19 patients compared to those with mild disease or healthy controls2. Recent studies, however, have found that patients with COVID-19 ARDS in fact have lower plasma cytokine levels compared to those with ARDS due to other causes3, highlighting a need to understand the underlying mechanisms of COVID-19 ARDS.
Clinical trials have demonstrated a significant mortality benefit for dexamethasone in COVID-19 patients with ARDS4, implicating a role for dysregulated inflammation in COVID-19 pathophysiology given the immunomodulatory effects of corticosteroids. In contrast, clinical trials of corticosteroids for ARDS prior to the SARS-CoV-2 pandemic have had mixed results, ranging from benefit to possible harm1. These differences suggest distinct biology in COVID-19 ARDS, with important implications for pathogenesis and treatment.
While several studies have assessed host airway transcriptional responses to SARS-CoV-25,6, none have compared COVID-19 ARDS to other causes of ARDS. We hypothesized that this comparison would identify distinct biological features of SARS-CoV-2 related lung injury, and thus conducted a case-control study of 52 critically ill patients requiring mechanical ventilation (Supplementary Table 1) for ARDS from COVID-19 (COVID-ARDS, n= 15), ARDS from other etiologies (Other-ARDS, n= 32), or for airway protection in the absence of pulmonary disease (No-ARDS, n = 5). Other ARDS etiologies included pneumonia, aspiration, sepsis, and transfusion reaction.
Patients were enrolled at two tertiary care hospitals in San Francisco, California under research protocols approved by the University of California San Francisco Institutional Review Board (Supplementary Methods). We excluded immunosuppressed patients to avoid confounding the measurement of host inflammatory responses (Supplementary Methods). Tracheal aspirate (TA) was collected within five days of intubation and underwent RNA sequencing (Supplementary Methods).
We compared TA gene expression between COVID-ARDS and Other-ARDS patients (Supplementary Methods, Figure 1a, Supplementary Data 2) and identified 696 differentially expressed genes at an adjusted P-value < 0.1, as well as differentially activated pathways using Ingenuity Pathway Analysis (IPA)7. Notably, we did not observe elevated expression of genes encoding canonical proinflammatory cytokines, such as IL-1 or IL-6, in COVID-ARDS compared to Other-ARDS. In fact, the IL-1, IL-6 and IL-8 signaling pathways were more highly activated in Other-ARDS, whereas COVID-ARDS patients had comparable inflammatory pathway activation to No-ARDS controls (Figure 1b, Supplementary Data 3). We also found attenuation of the proinflammatory HIF-1a and nitric oxide signaling pathways in COVID-ARDS compared to Other-ARDS patients. To relate these lower respiratory tract findings to systemic inflammatory responses, we also assessed circulating plasma cytokines. We found no difference in IL-6 or IL-8 levels, and lower concentrations of sTNFR1, in COVID-ARDS versus Other-ARDS patients (Figure 1c, Supplementary Data 4).
Evaluation of genes with the most significant expression differences in COVID-ARDS compared to Other-ARDS did, however, reveal several differences in genes regulating immunity and inflammation (Supplementary Data 2). For instance, among genes upregulated in COVID-ARDS, P2RY14 functions in purinergic receptor signaling to mediate inflammatory responses and its ligand UDP-glucose promotes neutrophil recruitment in the lung8. Conversely, ARG1, which promotes macrophage efferocytosis and inflammation resolution9, was downregulated in COVID-ARDS versus Other-ARDS patients.
At the pathway level, we observed activation of PTEN signaling in COVID-ARDS compared to both Other-ARDS and No-ARDS patients (Figure 1B, Supplementary Data 3). PTEN modulates both innate and adaptive immune responses by opposing the activity of PI3K10. Consistent with our observations, PTEN attenuates expression of certain cytokines while amplifying other innate immune responses in a manner that may promote injurious inflammation during respiratory infections11, suggesting a potentially pathologic role in COVID-ARDS. In silico prediction of cell type composition (Methods, Supplementary Figure 1, Supplementary Data 5) did not identify differences in lymphocyte, macrophage or neutrophil populations but did identify markedly decreased proportions of type-2 alveolar epithelial cells and increased proportions of goblet and ciliated cells in COVID-ARDS compared to Other-ARDS. This may reflect alveolar epithelial injury, airway remodeling, and/or preferential SARS-CoV-2 infection of cells with the highest expression of ACE2 and TMPRSS212.
To test the hypothesis that existing pharmaceuticals might counter the dysregulated transcriptional signature of COVID-19 related ARDS, we compared differentially expressed genes against the IPA database of 12,981 drug treatment-induced transcriptional signatures derived from human studies and cell culture experiments7 (Supplementary Methods, Figure 1c, Supplementary Data 6). Dexamethasone was the compound predicted to most significantly counter-regulate the genes expressed in COVID-ARDS patients compared to No-ARDS patients, and was also significant with respect to the Other-ARDS group. This finding was striking given that dexamethasone is the only drug thus far proven to confer a mortality benefit in patients with severe COVID-194. Granulocyte colony stimulating factor (G-CSF) was also significant, which is intriguing given that a recent clinical trial found a mortality benefit in COVID-19 patients treated with this agent13. Other corticosteroids (fluticasone, prednisolone) as well as omega-3 fatty acids (eicosapentenoic and docosahexaenoic acids) were also predicted to counteract the transcriptional profile of COVID-ARDS with respect to comparator groups and thus may represent possible therapeutic agents (Figure 1c, Supplementary Data 6).
As our analysis did not reveal evidence of a cytokine storm in COVID-19 ARDS, we hypothesized that other steroid-responsive pathways may be responsible for the therapeutic benefit of dexamethasone. Although commonly thought of as indiscriminate immunosuppressive medications, glucocorticoids affect diverse biological processes. We therefore proceeded to examine the genes comprising the transcriptional signature of COVID-ARDS that were also predicted to be regulated by dexamethasone (Supplementary Data 7). Interestingly, both dexamethasone and G-CSF were predicted to attenuate the expression of several genes highly upregulated in COVID-ARDS versus controls (e.g., P2YR14) as well as other genes with well-established roles in immunity, inflammation, and interferon responses. For instance, we found that COVID-ARDS patients had increased expression of the interferon-inducible and dexamethasone-regulated gene EPSTI1, which promotes M1 macrophage polarization14, and STAT1, which induces chemokine expression, regulates differentiation of hematopoietic cells, and promotes reactive oxygen species production15.
ARDS is a heterogeneous syndrome caused by diverse infectious and non-infectious insults1. To more precisely understand host response relationships between subtypes of ARDS, we performed secondary analyses comparing the transcriptional signature of COVID-ARDS without co-infections (n = 8) to that of ARDS caused exclusively by other viral (n = 4, Figure 2a, Supplementary Data 1) or bacterial (n = 9, Figure 2b) lower respiratory tract infections (LRTI) (Supplementary Data 8). COVID-ARDS was characterized by lower expression of proinflammatory signaling pathways compared to bacterial LRTI-associated ARDS (Figure 2c, Supplementary Data 9), but higher levels of the same pathways compared to viral LRTI-associated ARDS.
Although interferon-related gene expression was higher in COVID-ARDS compared to bacterial LRTIs and no-ARDS controls, it was markedly attenuated in ARDS patients with COVID-19 versus those with other viral LRTI (Figure 2d, Supplementary Data 10, Supplementary Methods). Prior studies found strong correlation between SARS-CoV-2 viral load and expression of interferon-stimulated genes (ISGs) in the upper respiratory tract of patients with mild disease, early during infection16. In contrast, in the lower respiratory tract of patients with severe disease, we observed decoupling of this relationship for several ISGs (Figures 2e-f, Supplementary Data 11), suggesting that a dysregulated interferon response in the lower respiratory tract may characterize severe COVID-19. This hypothesis is supported by recent findings of impaired interferon signaling in peripheral blood immune cells of patients with severe versus mild COVID-1917, and a recent report suggesting that a dysregulated interferon response may be a universal characteristic of severe viral infections18.
Discussion
Our results challenge the cytokine storm model of COVID-19 ARDS. Instead, we observe a complex picture of host immune dysregulation that includes upregulation of genes with non-canonical roles in inflammation, immunity and interferon signaling that are predicted to be attenuated by dexamethasone, G-CSF and other potential therapeutics. Our work emphasizes the value of including clinically relevant control patients in COVID-19 immunophenotyping studies and underscores that detection of elevated cytokine levels in the blood does not necessarily equate to causality in pathogenesis. Single cytokine blockade was attempted unsuccessfully in the past for treatment of sepsis19, which like COVID-19, is characterized by dysregulated host response to infection as well as significant biologically meaningful heterogeneity20.
This work builds on recent reports of dysregulated interferon responses in patients with severe COVID-19 pneumonia and suggests that decoupling of viral load from interferon signaling may be a relevant pathologic feature of severe disease. Further work in a larger cohort that also includes direct measurement of cytokine levels in the lower airway will be needed to validate these results. In conclusion, comparative lower respiratory transcriptional profiling of patients with COVID-19 and other ARDS etiologies did not find evidence of cytokine storm but did reveal potential therapeutic agents and novel gene targets that may mitigate the dysregulated host response of COVID-ARDS.
Data and Code Availability
Human gene counts for the samples generated in this study can be obtained under NCBI GEO accession GSE163426. The published human lung single-cell datasets17 used for cell type proportions analysis can be obtained through Synapse under accessions syn21560510 and syn21560511. Code for the differential expression and cell type proportions analysis is available at: https://github.com/AartikSarma/COVIDARDS.
Data Availability
Human gene counts for the samples generated in this study can be obtained under NCBI GEO accession GSE163426. The published human lung single-cell datasets used for cell type proportions analysis can be obtained through Synapse under accessions syn21560510 and syn21560511. Code for the differential expression and cell type proportions analysis is available on GitHub (see link).
Author Contributions
CRL, CSC, AS and SC conceived and designed the study. TD, RG, PHS and KMA oversaw or performed sample processing, library preparation and sequencing. CSC, CMH, KNK RG, AJ, JGW and ERS coordinated or contributed to clinical operations including patient enrollment. CD, AS, CRL, FM, TD and NS performed metadata collation or clinical chart review. AS, EM, AOP, CRL and SAC performed data analysis and interpretation. CSC, SAC, EM, DJE, JLD, KMA, MFK, PFW, MAM, BSZ, JGW, AL, AB, FM, PS and MS provided guidance, advice and comments on the study design and manuscript. CRL, AS and CSC wrote the manuscript with input from all authors.
Supplementary Materials
Supplementary Methods
Study design, clinical cohort and ethics statement
We conducted a case-control study of patients with ARDS due to COVID-19 (n = 15) versus two control groups of either patients with ARDS due to other causes (n = 32) or patients intubated for airway protection without evidence of pulmonary pathology on imaging (n = 5). We studied patients who were enrolled in either of two prospective cohort studies of critically ill patients at the University of California, San Francisco (UCSF) and Zuckerberg San Francisco General Hospital. Both studies were approved by the UCSF Institutional Review Board according under protocols 17-24056 and 20-30497, respectively, which granted a waiver of initial consent for tracheal aspirate and blood sampling. Informed consent was subsequently obtained from patients or their surrogates for continued study participation, as previously described1.
For this analysis, inclusion criteria were: 1) admission to the intensive care unit for mechanical ventilation for ARDS or airway protection, 2) age ≥ 18 years, 3) availability of TA with 106 protein-coding reads collected within five days of intubation. Exclusion criteria were: 1) receipt of immunosuppressive medication or underlying immunocompromising condition prior to tracheal aspirate collection. Subject charts and chest x-rays were reviewed by at least two study authors (AS, PS, ES, FM, CD, MM, CL, CC) to confirm a diagnosis of ARDS using the Berlin Definition2. Lower respiratory tract infections were adjudicated by two study physicians using the United States Centers for Disease Control surveillance definition of pneumonia3. Of 72 patients initially screened, nine with ARDS due to COVID-19 and 10 with ARDS due to other causes were excluded because of treatment with immunosuppressant medications or because of an underlying immunocompromising condition (e.g., solid organ transplantation, bone marrow transplantation, human immunodeficiency virus infection).
Metagenomic sequencing
Following enrollment, TA was collected and mixed 1:1 with DNA/RNA shield (Zymo Research) to preserve nucleic acid. To evaluate host gene expression and detect the presence of SARS-CoV-2 and other viruses, metagenomic next generation sequencing (mNGS) of RNA was performed on TA specimens. Following RNA extraction (Zymo Pathogen Magbead Kit) and DNase treatment, human cytosolic and mitochondrial ribosomal RNA was depleted using FastSelect (Qiagen). To control for background contamination, we included negative controls (water and HeLa cell RNA) as well as positive controls (spike-in RNA standards from the External RNA Controls Consortium (ERCC))4. RNA was then fragmented and underwent library preparation using the NEBNext Ultra II RNAseq Kit (New England Biolabs). Libraries underwent 146 nucleotide paired-end Illumina sequencing on an Illumina Novaseq 6000 instrument.
Host differential expression and pathway analysis
Following demultiplexing, sequencing reads were pseudo-aligned with kallisto5 (v. 0.46.1; including bias correction) to an index consisting of all transcripts associated with human protein coding genes (ENSEMBL v. 99), cytosolic and mitochondrial ribosomal RNA sequences, and the sequences of ERCC RNA standards. Samples retained in the dataset had a total of at least 1,000,000 estimated counts associated with transcripts of protein coding genes, and the median across all samples was 7.3 million. Gene-level counts were generated from the transcript-level abundance estimates using the R package tximport6, with the scaledTPM method.
Differential expression analysis was performed using DESeq27. We modeled the expression of individual genes using the design formula ∼ARDSEtiology. In our primary analysis, the ARDS etiology was categorized as COVID-ARDS, Other-ARDS, or No-ARDS. In our secondary analysis, the ARDS etiology was categorized as COVID-ARDS, Viral-ARDS, Bacterial-ARDS, or No-ARDS. COVID-ARDS patients with viral or bacterial co-infections were excluded from this secondary analysis. Significant genes were identified using an independent-hypothesis-weighted, Benjamini-Hochberg false discovery rate (FDR) less than 0.18,9. Empirical Bayesian adaptive shrinkage estimators for log2-fold change were fit using ashr10. We generated heatmaps of the top 50 differentially expressed genes by absolute log2-fold change. For visualization, gene expression was normalized using the variance stabilizing transformation, centered, and z-scaled. Heatmaps were generated using the pheatmap package. Patients were clustered using Euclidean distance and genes were clustered using Manhattan distance. Differentially expressed genes (FDR < 0.1) were analyzed using Ingenuity Pathway Analysis (IPA, Qiagen)11,12.
Canonical pathway analysis and drug/cytokine upstream regulator analysis
To evaluate signaling pathways and upstream transcriptional regulators from gene expression data, we employed IPA. Specifically, genes were analyzed using Core, Canonical Pathway and Upstream Regulator Analysis on shrunken log2-fold change. IPA Upstream Regulator Analysis was employed to identify potential drug and cytokine regulators and predict their activation states based on expected effects between regulators and their known target genes or proteins annotated in the Ingenuity Knowledge Base (IKB)12. IPA calculates a Fisher’s exact p-value for overlap of differentially expressed genes with curated gene sets representing canonical biological pathways, or upstream regulators of gene expression, including cytokines and 12,981 drugs. In addition, IPA calculates a Z-score for the direction of gene expression for a pathway or regulator based on the observed gene expression in the dataset. The Z-score signifies whether expression changes for genes within pathways, or for known target genes of upstream regulators, are consistent with what is expected based on previously published analyses annotated in the IKB. Significant pathways and upstream regulators were defined as those with a Z-score absolute value greater than 2 and an overlap P value < 0.05.
In silico analysis of cell type proportions
Cell-type proportions were estimated from bulk host transcriptome data using the CIBERSORT X algorithm13. We used the Human Lung Cell Atlas dataset14 to derive the single cell signatures. The cell types estimated with this reference cover all expected cell types in the airway. The estimated proportions were compared between the three patient groups using a Mann-Whitney-Wilcoxon test (two-sided) with Bonferroni correction.
Quantification of SARS-CoV-2 viral load by mNGS
All samples were processed through a SARS-CoV-2 reference-based assembly pipeline that involved removing reads likely originating from the human genome or from other viral genomes annotated in RefSeq with Kraken2 (v. 2.0.8_beta), and then aligning the remaining reads to the SARS-CoV-2 reference genome MN908947.3 using minimap2 (v. 2.17). We calculated SARS-CoV-2 reads-per-million (rpM) by dividing the number of reads that aligned to the virus with mapq ≥ 20 by the total number of reads in the sample (excluding reads mapping to ERCC RNA standards).
Regression of ISG counts against viral load in TA and NP samples
We assembled a set of 100 interferon-stimulated genes based on the “Hallmark interferon alpha response” gene set in MSigDB15. We then performed robust regression of the quantile normalized gene counts (log2 scale), generated using the R package limma, against log10(rpM) of SARS-CoV-2. This was done in two separate datasets of COVID-19 patients: i) the tracheal aspirate (TA) samples from patients with COVID-19 ARDS reported in this study (n=15); and ii) the nasopharyngeal swab (NP) samples from patients with mostly early and mild disease that we previously reported (n=93)16. The analysis was performed using the R package robustbase (v. 0.93.6), which implements MM-type estimators for linear regression. Model predictions were generated using the R package ggeffects (v. 0.14.3) and used for display in the individual gene plots. Error bands represent normal distribution 95% confidence intervals around each prediction. Reported P-values for significance of the difference of the regression coefficient from 0 are based on a t-statistic and Benjamini-Hochberg adjusted. Reported R2 values represent the adjusted robust coefficient of determination.
Data and Code Availability
Human gene counts for the samples generated in this study can be obtained under NCBI GEO accession GSE163426. The published human lung single-cell datasets17 used for cell type proportions analysis can be obtained through Synapse under accessions syn21560510 and syn21560511. Code for the differential expression analyses, cell type proportions analysis and gene expression classifiers is available at: https://github.com/AartikSarma/COVIDARDS.
Supplementary Data 1. Detailed clinical, microbiological and demographic features of the cohort. Legend: *SaO2/FiO2 Ratio (SF Ratio) < 315 was used to verify ARDS diagnosis in these subjects.
Supplementary Data 2. Differentially expressed genes (adjusted P value (padj) < 0.1) between patients with ARDS due to COVID-19 (COVID-ARDS) versus a) controls with ARDS due to other etiologies (Other-ARDS) or b) intubated controls without ARDS (No-ARDS). Positive fold change indicates gene is upregulated in COVID-ARDS.
Supplementary Data 3. Ingenuity Pathway Analysis (IPA) of differentially expressed genes (padj < 0.1) between patients with ARDS due to COVID-19 (COVID-ARDS) versus a) controls with ARDS due to other etiologies (Other-ARDS) or b) intubated controls without ARDS (No-ARDS). Positive Z-score indicates pathway is upregulated in COVID-ARDS. Pathways with a Z-score absolute value ≥ 1 are included in table.
Supplementary Data 4. Plasma concentrations of inflammatory cytokines in patients with ARDS due to COVID-19 (COVID-ARDS) versus controls with ARDS due to other etiologies (Other-ARDS).
Supplementary Data 5. In silico deconvolution of cell type proportions from tracheal aspirate bulk RNA-sequencing data using lung single cell signatures. Data are plotted in (Supplementary Figure 1).
Supplementary Data 6. Chemical and biological drugs computationally predicted by IPA to attenuate the transcriptional response of ARDS due to COVID-19 (COVID-ARDS) against a comparator group of a) ARDS due to other causes (Other-ARDS) or b) intubated controls without ARDS (No-ARDS). Drugs with a Z-score > 2 are included in table.
Supplementary Data 7. Genes affected by drugs computationally predicted by IPA to modulate the transcriptional response of ARDS due to COVID-19 (COVID-ARDS) against a comparator group of a) ARDS due to other causes (Other-ARDS) or b) intubated controls without ARDS (No-ARDS). Predicted transcriptional effect of each drug on each gene is indicated with respect to comparator group 1 in the table.
Supplementary Data 8. Differentially expressed genes (padj < 0.1) between patients with ARDS due to COVID-19 (COVID-ARDS) versus controls with ARDS due to a) other viral lower respiratory tract infections (Other Viral-ARDS) or b) bacterial lower respiratory tract infections (Bacterial-ARDS). Positive log2 fold change indicates gene is upregulated in COVID-ARDS.
Supplementary Data 9. Pathway analysis (IPA) of differentially expressed genes (padj < 0.1) between patients with ARDS due to COVID-19 (COVID-ARDS) versus controls with ARDS due to other viral or bacterial lower respiratory tract infections. Z-scores are with respect to COVID-ARDS. Pathways with a Z-score absolute value ≥ 1 are included in table.
Supplementary Data 10. Computationally predicted (IPA) upstream cytokines based on the transcriptional signature of COVID-19 ARDS (COVID-ARDS) compared to ARDS from other viral LRTI. Pathways with a Z-score absolute value ≥ 1 are included in table.
Supplementary Data 11. Relationship between interferon-stimulated gene (ISG) expression and SARS-CoV-2 viral load measured in RNA-seq reads per million (rpM) for COVID-19 patients with late, severe disease and ARDS (COVID-ARDS) from lower respiratory tract samples (this study), and COVID-19 patients with early, mostly mild disease measured from upper respiratory tract samples16. Legend: reg_intercept = intercept in the robust regression of gene expression on SARS-CoV-2 viral load; reg_slope = slope in the robust regression of gene expression on SARS-CoV-2 viral load; reg_adj_R2 = adjusted robust coefficient of determination; reg_p_adj = Benjamini-Hochberg adjusted p-value for difference of the regression slope from 0.
Supplementary Appendix. COMET Consortium Member list.
Acknowledgements
This study was performed with support from the National Institute of Allergy and Infectious Diseases-sponsored Immunophenotyping Assessment in a COVID-19 Cohort (IMPACC) Network. We gratefully appreciate support from Amy Kistler, PhD, Jack Kamm, PhD, Saharai Caldera, BS, Maira Phelps BS and Norma Neff, PhD. We thank Mark and Carrie Casey, Julia and Kevin Hartz, Carl Kawaja and Wendy Holcombe, Eric Keisman and Linda Nevin, Martin and Leesa Romo, Three Sisters Foundation, Diana Wagner and Jerry Yang and Akiko Yamazaki for their support.
Footnotes
↵+ COMET (COVID-19 Multiphenotyping for Effective Therapies) Consortium members are listed in the Supplementary Appendix.
Funding NIH R35 HL140026 (CSC), NIH/NIAID U19 AI1077439 (DE), NHLBI K23HL138461-01A1 (CL), K24HL137013 (PGW), F32 HL151117 (AS), Chan Zuckerberg Biohub (AOP, AB, JLD), Philanthropic contributions from: Mark and Carrie Casey, Julia and Kevin Hartz, Carl Kawaja and Wendy Holcombe, Eric Keisman and Linda Nevin, Martin and Leesa Romo, Three Sisters Foundation, Diana Wagner, Jerry Yang and Akiko Yamazaki.
References
Supplementary References
- 1.
- 2.
- 3.
- 4.
- 5.
- 6.
- 7.
- 8.
- 9.
- 10.
- 11.
- 12.
- 13.
- 14.
- 15.
- 16.
- 17.