key: cord-0873564-ehgurh2w authors: Overholt, Kalon J.; Krog, Jonathan R.; Zanoni, Ivan; Bryson, Bryan D. title: Dissecting the common and compartment-specific features of COVID-19 severity in the lung and periphery with single-cell resolution date: 2021-06-17 journal: iScience DOI: 10.1016/j.isci.2021.102738 sha: d40bb56c12beed1bb64cc858f6b7a0a4972844c3 doc_id: 873564 cord_uid: ehgurh2w Severe COVID-19 is accompanied by rampant immune dysregulation in the lung and periphery, with immune cells of both compartments contributing to systemic distress. The extent to which immune cells of the lung and blood enter similar or distinct pathological states during severe disease remains unknown. Here, we leveraged 96 publicly available single-cell RNA sequencing datasets to elucidate common and compartment-specific features of severe-to-critical COVID-19 at the levels of transcript expression, biological pathways, and ligand-receptor signaling networks. Comparing severe patients to milder and healthy donors, we identified distinct differential gene expression signatures between compartments but a core set of co-directionally regulated surface markers. A majority of severity-enriched pathways were shared, while TNF and interferon responses were polarized. Severity-specific ligand-receptor networks appeared to be differentially active in both compartments. Overall, our results describe a nuanced response during severe COVID-19 where compartment plays a role in dictating the pathological state of immune cells. patients showed only a few commonalities between the lung and periphery. This generally 158 dissimilar trend followed when using an alternative approach to measure dissimilarity by directly 159 conducting differential gene expression between cells of different compartments ( Figure S4A , 160 see STAR Methods). 161 To identify the most strongly evidenced differentially expressed transcripts common to 162 both compartments, we used a DEG classification scheme with more stringent p-value and fold 163 change thresholds and explicitly accounted for donor-to-donor variation in the statistical test 164 (see STAR Methods). The stringently selected DEGs that were observed in both the Liao BALF 165 dataset and at least 2 PBMC datasets are given in Supplemental Table S1 . Declines in HLA 166 class II related gene expression (HLA-DPs, HLA-DRs), upregulation of calprotectin and 167 calgranulin genes (S100A8 and S10012), and upregulation of the clusterin gene CLU were 168 common features of severe patients in both the lung and periphery. pathways, indicating that these pathways are generally active in severe disease across nearly 208 every immune cell type ( Figure 4A ). In the Lee PBMC cohort, many of the same enriched 209 pathways were observed, including IFN-α response, IFN-γ response, and TNF-α signaling via 210 NF-κB ( Figure 4A) . Interestingly, the majority of differentially regulated pathways appeared to 211 be conserved between the lung and blood for most cell types (Figure 4B) , and this result was 212 consistent across cohorts ( Figure 4B and Figure S7 ). 213 To examine pathway-level regulation specific to severe disease courses, we applied the 214 same analysis to DEGs between severe and mild COVID-19 patients ( Figure 4C ). In the lung, a 215 number of broadly enriched pathways including IL-2/STAT5 signaling and TNF-α signaling via 216 NF-kB were observed ( Figure 4C) . Strikingly, type I/III and type II interferon responses showed 217 mixed enrichment and depletion across lung cell types, with enrichment tending to occur in 218 myeloid cells (CD14 and CD16 monocytes, cDC2s, neutrophils) and depletion tending to occur 219 in lymphoid cells (T cells, B cells, and pDCs). In the blood, more consistent enrichment of IFN-α 220 J o u r n a l P r e -p r o o f and IFN-γ responses was observed, without a clear distinction between myeloid cells and 221 lymphoid cells ( Figure 4C) . In contrast to the lung, the TNF-α signaling pathway was 222 significantly depleted in multiple blood cell types including CD4 Treg, CD8 effector and naïve T, 223 and NK cells. Pooling data from all the available cohorts, we observed a similar effect in which 224 the TNF pathway was more strongly active in the blood compared to lung for healthy and mild 225 donors, but more active in the lung than in the blood for severe donors (Figure S4B , see STAR 226 methods). As many of the cell types studied via GSEA showed significant enrichment of both 227 the "interferon alpha response" and "interferon gamma response" hallmark pathways, we sought 228 to identify whether these pathways had been induced by the same gene sets. The degree of 229 overlap (Jaccard index of the GSEA leading edge, see STAR Methods) was less than 50% in 230 the large majority of cell types and conditions ( Figure S8 ). On the whole, alterations in pathway-231 level activity differentiating severe from mild patients demonstrated greater than 20% 232 conservation in most cell types across the cohorts compared in this analysis, but a substantial 233 fraction of contra-directionally regulated or "polarized" pathways including TNF-α signaling and 234 interferon responses were detected across all of the cohorts studied ( Figure 4D and Figure information to predict whether transcriptional signatures were likely induced through intra-245 compartmental signaling or whether the data were consistent with cross-compartmental 246 signaling, termed "crosstalk" (see Discussion for alternative mechanisms). Briefly, differentially 247 expressed ligands (DELs) were identified using NicheNet through: 1) a data-driven unbiased 248 approach searching for ligands differentially expressed between the conditions of interest that 249 act broadly to induce differential gene expression in over one-third of cell types in a receiver 250 compartment, and 2) a targeted approach to find specific transcripts of interest (see STAR 251 Methods). Following the identification of DELs, we categorized each according to whether their 252 "footprint" in a given compartment could be explained by the DEL's upregulation in the same 253 compartment (suggesting intra-compartmental signaling), the opposite compartment 254 (suggesting crosstalk) or both compartments. The DEL profile for differential gene expression between severe patients and healthy 256 controls is shown in Figure 5A and is schematically represented in Figure 5B . TNF-induced 257 differential expression signatures along with upregulation of TNF occurred in both 258 compartments, suggesting intra-compartmental signaling. Upregulation of genes including 259 TGFB1, CCL2, and SPP1 was observed only in the lung. However, over one-third of cell types 260 in the blood bore differential gene expression signatures induced by these ligands, indicating 261 potential lung-to-blood crosstalk via these secreted factors. On the other hand, gene expression 262 signatures downstream of IFN-γ and IL-1α were widespread in the lung, but IFNG and IL1A 263 were only upregulated in the blood, suggesting blood-to-lung crosstalk. Using the same analysis pipeline to examine differential gene expression programs 265 between severe and mild disease ( Figure 5C , schematically represented in Figure 5D ) showed We identified a small but persistent degree of overlap between comparable cell types in The study is subject to several limitations. First, our study does not include a cohort in 352 which cells of the lung and the blood have been analyzed in the same patients although this 353 type of cohort would be ideal for validating our predictions. Next, the ability to use healthy BALF 354 data as a reference point remains limited without conducting differential gene expression tests 355 between cohorts collected by separate investigators in separate studies. Additionally, the set 356 publicly available mild COVID-19 patient BALF scRNA-seq samples remains small and limits 357 our ability to make entirely orthogonal comparisons between the lung and blood in this study. Finally, this study does not contain experimental validation of the putative ligand-receptor 359 networks predicted to be active in the lung and blood or the potential crosstalk between these 360 networks. Figure S6 and Table S1 . Gene-barcode matrices obtained from the GEO were preprocessed using Seurat (v. Re-integration was peformed using Harmony as explained above, but for myeloid cells the PCs 572 were calculated using the point where the first difference in standard deviation was less than 573 0.1% to produce the cleanest clustering result. Cell types were then annotated using scaled 574 expression levels of canonical marker genes as explained above. Figure S1 shows the levels of 575 key marker genes for the sub-cluster analysis. We took multiple steps to improve the purity of sub-clusters. Following cell type 579 annotation, clusters labeled as putative doublets (positive for multiple lineage markers) were 580 removed from the analysis. Additionally, cells containing more than 5% total hemoglobin (HBA, 581 HBB, HBD) transcripts or more than 1% PPBP transcripts were considered to be either (Table S1 ) 606 To identify the highest-confidence DEGs using stringent methods, the same DE analysis Figure S4) . 616 In the supplementary analyses of Figure S4 , we desired to study transcriptional 617 differences between compartments for a given cell type in a given severity condition (e.g. severe neutrophils in BALF vs. severe neutrophils in PBMCs). To do this, it was necessary to 619 conduct differential gene expression tests between donors who had been sequenced as part of 620 different studies. We used the following method to attempt to mitigate the issue of batch effect. First, all BALF donors from the severity condition of interest were combined into a single pool 622 and all PBMC donors from that condition were combined into a separate pool. Next, we 623 conducted differential gene expression between these pools in a given cell type using the 624 'FindMarkers' function in Seurat implementing the MAST algorithm. In addition to using the UMI 625 count depth as a numerical covariate in MAST, we also used the study of origin as a categorical 626 covariate in the 'latent.vars' variable of the 'FindMarkers' function. This allowed us to explicitly 627 model the batch effects between studies using a fixed effects model. As a result, the DEG list 628 was smaller than using the same approach without considering the study of origin as a 629 covariate. In Figure S4 , we considered DEGs that were robust to study-to-study batch variation leading edge subsets for both pathways were compared to identify unique and common genes. Overlap was quantified using the Jaccard similarity index (intersection of leading edge subsets/union of leading edge subsets). The Jaccard index and the genes uniquely contributing 658 to either pathway were visualized for each cell type ( Figure S8 ). Compartmental Comparisons for Main Pipeline (Figures 2-5 ) 661 In Figures 2-5 or significantly downregulated in both compartments was identified and termed "co-directional". The intersecting set of genes or pathways that were significantly regulated in different directions 668 between compartments was identified and termed "contradirectional". The remaining genes not 669 part of the intersection were termed "BALF only" or "PBMC only". A modified Jaccard similarity To generate the rose plots in Figures 2 and 3 , we considered only DEGs with adjusted However, we did not filter these genes prior to conducting GSEA as the GSEA algorithm 684 requires the total gene set as input (Subramanian et al., 2005) . Analysis of Between-Compartment DEGs ( Figure S4 ) 687 In Figure S4 , we established DEGs between compartments for a given cell type in a 688 given severity condition, giving for example a differential expressed gene list between severe 689 neutrophils in the lung and severe neutrophils in blood. These DEG lists were first filtered as TGFB1 TNF TNFSF14 IFNG IL1B TNF TNFSF10 TNFSF14 TNFSF18 ADAM17 ADAM9 CD99 CFH CTGF GZMB ICAM1 ICAM3 IFNG IL15 ITGB2 LAMC2 SPP1 TGFB1 THBS1 TNF IFNG TNF TNFSF10 TNFSF13B TNFSF14 TNFSF8 ADAM17 CCL2 CCL3 CCL8 IFNG IL15 SPP1 TGFB1 TNF IFNG TNF TNFSF10 TNFSF13B IFNG IL15 CCL2 CCL3 CCL8 IFNG IL15 SPP1 TGFB1 TNFSF14 TNF TNF CCL2 CCL8 IL15 SPP1 TGFB1 LIF TGFB1 TNFSF14 IFNG IL1A IL15 SPP1 TGFB1  96 single-cell public datasets from severe, mild-to-moderate, and healthy donors  Severity-specific cell surface markers and pathways shared between compartments  Putative ligand-receptor signaling dialogs within and between compartments Imbalanced Host Response to SARS-CoV-2 Drives 801 Development of COVID-19 Host-Viral Infection Maps Reveal Signatures of Severe COVID-19 Patients Type III interferons disrupt the 809 lung epithelial barrier upon viral recognition NicheNet: modeling intercellular communication 812 by linking ligands to target genes A familial cluster of 817 pneumonia associated with the 2019 novel coronavirus indicating person-to-person 818 transmission: a study of a family cluster Transcriptomic Analysis of COVID-19 Blood, Lung, and Airway. BioRxiv, Microbe Interferon responses in viral pneumonias Baseline 851 Characteristics and Outcomes of 1591 Patients Infected with SARS-CoV-2 Admitted Impaired 857 type I interferon activity and inflammatory responses in severe COVID-19 patients Clinical features of patients infected with 2019 novel coronavirus in Wuhan, China. The 862 Cross-country comparison of case fatality rates of Covid-864 Fast gene set enrichment analysis Gene 941 set enrichment analysis: A knowledge-based approach for interpreting genome-wide 942 expression profiles A molecular cell atlas of the 947 human lung from single-cell RNA sequencing Discriminating mild from critical COVID-19 by innate and 953 adaptive immune single-cell profiling of bronchoalveolar lavages A single-cell atlas of the peripheral immune response severe COVID-19 Cytokine profile in plasma of severe COVID-963 19 does not differ from ARDS and sepsis A new coronavirus associated with human 968 respiratory disease in China Characteristics of and Important Lessons from the 971 COVID-19) Outbreak in China: Summary of a Report of 72314 Transcriptomic characteristics of bronchoalveolar lavage fluid and peripheral blood 978 mononuclear cells in COVID-19 patients The differential immune responses to COVID-19 in 982 peripheral and lung revealed by single-cell RNA sequencing Pathological findings of COVID-19 associated with acute respiratory distress syndrome Clinical 991 course and risk factors for mortality of adult inpatients with COVID-19 in Wuhan, China: a 992 retrospective cohort study The authors gratefully acknowledge Paul Blainey, Cal Gunnarsson, Bianca Lepe, Megan Tse, 364