key: cord-0989284-x6akbo5f
authors: Sarma, A.; Christenson, S. A.; Mick, E.; DeVoe, C.; Deiss, T.; Pisco, A. O.; Ghale, R.; Jauregui, A.; Byrne, A.; Moazed, F.; Spottiswoode, N.; Sinha, P.; Zha, B. S.; Serpa, P. H.; Ansel, K. M.; Wilson, J. G.; Leligdowicz, A.; Siegel, E. R.; Sirota, M.; DeRisi, J. L.; Matthay, M. A.; COMET Consortium,; Hendrickson, C. M.; Kangelaris, K. N.; Krummel, M.; Woodruff, P. G.; Erle, D. J.; Calfee, C. S.; Langelier, C. R.
title: COVID-19 ARDS is characterized by a dysregulated host response that differs from cytokine storm and may be modified by dexamethasone
date: 2021-01-04
journal: nan
DOI: 10.1101/2020.12.28.20248552
sha: f989af7d560450a128c35dd04dadb1303838d2d2
doc_id: 989284
cord_uid: x6akbo5f

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.

with non-canonical roles in inflammation and immunity that were predicted to be modulated by CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted January 4, 2021. ; https://doi.org/10.1101/2020.12.28.20248552 doi: medRxiv preprint

In its most severe form, coronavirus disease 2019 (COVID-19) can precipitate the acute 56 respiratory distress syndrome (ARDS), which is characterized by low arterial oxygen 57 concentrations, alveolar injury and a dysregulated inflammatory response in the lungs 1 CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted January 4, 2021. Ingenuity Pathway Analysis (IPA) 7 . Notably, we did not observe elevated expression of genes 88 encoding canonical proinflammatory cytokines, such as IL-1 or IL-6, in COVID-ARDS compared 89 to Other-ARDS. In fact, the IL-1, IL-6 and IL-8 signaling pathways were more highly activated in 90 Other-ARDS, whereas COVID-ARDS patients had comparable inflammatory pathway activation 91 to No-ARDS controls (Figure 1b, Supplementary Data 3) . We also found attenuation of the 92 proinflammatory HIF-1a and nitric oxide signaling pathways in COVID-ARDS compared to Other-93 ARDS patients. To relate these lower respiratory tract findings to systemic inflammatory 94 responses, we also assessed circulating plasma cytokines. We found no difference in IL-6 or IL- 

Evaluation of genes with the most significant expression differences in COVID-ARDS 98 compared to Other-ARDS did, however, reveal several differences in genes regulating immunity 

. CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review)

The copyright holder for this preprint this version posted January 4, 2021. 

and was also significant with respect to the Other-ARDS group. This finding was striking given 123 that dexamethasone is the only drug thus far proven to confer a mortality benefit in patients with 124 severe COVID-19 4 . Granulocyte colony stimulating factor (G-CSF) was also significant, which is 125 intriguing given that a recent clinical trial found a mortality benefit in COVID-19 patients treated 126 with this agent 13 . Other corticosteroids (fluticasone, prednisolone) as well as omega-3 fatty acids 127 (eicosapentenoic and docosahexaenoic acids) were also predicted to counteract the 

As our analysis did not reveal evidence of a cytokine storm in COVID-19 ARDS, we 131 hypothesized that other steroid-responsive pathways may be responsible for the therapeutic 132 . CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review)

The copyright holder for this preprint this version posted January 4, 2021. ; https://doi.org/10.1101/2020.12.28.20248552 doi: medRxiv preprint benefit of dexamethasone. Although commonly thought of as indiscriminate immunosuppressive 133 medications, glucocorticoids affect diverse biological processes. We therefore proceeded to 134 examine the genes comprising the transcriptional signature of COVID-ARDS that were also 135 predicted to be regulated by dexamethasone (Supplementary Data 7) . Interestingly, both 136 dexamethasone and G-CSF were predicted to attenuate the expression of several genes highly 137 upregulated in COVID-ARDS versus controls (e.g., P2YR14) as well as other genes with well- 

. CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted January 4, 2021. ; https://doi.org/10.1101/2020.12.28.20248552 doi: medRxiv preprint 209 210 

a b d PRKN FBN2 MZB1 GABRB3 DEFB115 LGI1 FRK OR4F17 TBC1D8B ADAMTSL3 NMNAT3 PIGL AMDHD1 LRRTM4 AC005726.1 SLC15A1 PTRH1 SLC12A5 SSBP2 WDR72 NEDD4 ATP1A3 DUOXA2 DUOX2 GSDMB MPRIP MX1 HERC6 USP18 CSF3 PML ISG15 PRKD2 RUBCN HELZ2 AC118281.1 ZNF213 IRGQ COLEC12 KRTAP4−16 ADAD2 MOBP MYCBP2 BAG1 F8A2 NTNG2 ZNF264 ZNF580 TNK2 CSNK1E DDX3Y ZNF354B TNFSF15 AC022506.1 GBE1 SGTB FCHSD2 USF3 TRIM33 ZFYVE16 UBE2B ZNF438 RIT1 STK4 RNF13 SDCBP ST20 AFF2 RIMBP3C STEAP4 CAPZA2 QPCT CACUL1 RNF141 CASP5 KRAS C9orf72 RAB32 ASAH1 RNF130 TMSB4X TGFBR1 LAMP2 ATP6V1C1 ABHD5 YPEL5 SGK1 CXCR4 CREB5 LTA4H APC TNFSF8 LSMEM1 AHCTF1 USP9X OSER1 MTHFS H3−3A MAP1LC3B CLEC7A dx

Other-viral Non-viral logFC -5 0 5

Interferon-λ1

Interferon-α1

Interferon-β1

Interferon-α2 IPA Z-score Increased Decreased COVID-ARDS vs. 

. CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted January 4, 2021. is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted January 4, 2021. ; https://doi.org/10.1101/2020.12.28.20248552 doi: medRxiv preprint

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 described 1 .

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 10 6 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 Definition 2 . Lower respiratory tract infections were adjudicated by two study physicians using the United States Centers for Disease Control surveillance definition of pneumonia 3 . 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 . CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted January 4, 2021. ; https://doi.org/10.1101/2020.12.28.20248552 doi: medRxiv preprint underlying immunocompromising condition (e.g., solid organ transplantation, bone marrow transplantation, human immunodeficiency virus infection).

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 kallisto 5 (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 tximport 6 , with the scaledTPM method.

Differential expression analysis was performed using DESeq2 7 . We modeled the expression of individual genes using the design formula ~ARDSEtiology. In our primary . CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted January 4, 2021. ; https://doi.org/10.1101/2020.12.28.20248552 doi: medRxiv preprint 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 independenthypothesis-weighted, Benjamini-Hochberg false discovery rate (FDR) less than 0.1 8, 9 . Empirical Bayesian adaptive shrinkage estimators for log2-fold change were fit using ashr 10 . 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 is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted January 4, 2021. ; https://doi.org/10.1101/2020.12.28.20248552 doi: medRxiv preprint 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 algorithm 13 . We used the Human Lung Cell Atlas dataset 14 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).

We assembled a set of 100 interferon-stimulated genes based on the "Hallmark interferon alpha response" gene set in MSigDB 15 . 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 . CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted January 4, 2021. ; https://doi.org/10.1101/2020.12.28.20248552 doi: medRxiv preprint 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 R 2 values represent the adjusted robust coefficient of determination.

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 17 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. . CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review)

The copyright holder for this preprint this version posted January 4, 2021. . CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review)

The copyright holder for this preprint this version posted January 4, 2021. ; https://doi.org/10.1101/2020.12.28.20248552 doi: medRxiv preprint . CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review)

The copyright holder for this preprint this version posted January 4, 2021. ; https://doi.org/10.1101/2020.12.28.20248552 doi: medRxiv preprint

Integrating host response and unbiased microbe detection for lower respiratory tract infection diagnosis in critically ill adults

Acute respiratory distress syndrome: the Berlin definition

United States Centers for Disease Control and Prevention. CDC/NHSN Surveillance Definitions for Specific Types of Infections

Evaluation of the External RNA Controls Consortium (ERCC) reference material using a modified Latin square design

Near-optimal probabilistic RNA-seq quantification

Differential analyses for RNA-seq: transcriptlevel estimates improve gene-level inferences

Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2

Data-driven hypothesis weighting increases detection power in genome-scale multiple testing

Controlling the False Discovery Rate: A Practical and Powerful Approach to Multiple Testing

False discovery rates: a new deal

Causal analysis approaches in Ingenuity Pathway Analysis

International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review) The copyright holder for this preprint this version posted

Robust enumeration of cell subsets from tissue expression profiles

A molecular cell atlas of the human lung from single-cell RNA sequencing

Gene set enrichment analysis: A knowledge-based approach for interpreting genome-wide expression profiles

Upper airway gene expression reveals suppressed immune responses to SARS-CoV-2 compared with other respiratory viruses

A molecular cell atlas of the human lung from single cell RNA sequencing

International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity

In silico deconvolution of cell types from tracheal aspirate bulk RNAsequencing data using lung single cell signatures. The horizontal line inside the box denotes the median and the lower and upper hinges correspond to the first and third quartiles, respectively.Whiskers extend from the hinge to the largest (smallest, respectively) value no more than 1.5*IQR away from the hinge, where IQR is the interquartile range. The y-axis in each panel was trimmed at the maximum value among the three patient groups of 1.5*IQR above the third quartile. Pairwise comparisons between patient groups were performed with a two-sided Mann-Whitney-Wilcoxon test followed by Bonferroni's correction (n=15 COVID-ARDS, n=32 Other ARDS, n=5 No-ARDS). Data are tabulated in (Supplementary Data 5).. CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review)The copyright holder for this preprint this version posted January 4, 2021. is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted January 4, 2021. ; https://doi.org/10.1101/2020.12.28.20248552 doi: medRxiv preprint samples (this study), and COVID-19 patients with early, mostly mild disease measured from upper respiratory tract samples 16 . CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review)The copyright holder for this preprint this version posted January 4, 2021. ; https://doi.org/10.1101/2020.12.28.20248552 doi: medRxiv preprintThe copyright holder for this preprint this version posted January 4, 2021. ; https://doi.org/10.1101/2020.12.28.20248552 doi: medRxiv preprint