Cellular and molecular heterogeneities and signatures, and pathological trajectories of fatal COVID-19 lungs defined by spatial single-cell transcriptome analysis ================================================================================================================================================================== * Arun Das * Wen Meng * Zhentao Liu * Md Musaddaqul Hasib * Hugh Galloway * Suzane Ramos da Silva * Luping Chen * Gabriel L Sica * Alberto Paniz-Mondolfi * Clare Bryce * Zachary Grimes * Emilia Mia Sordillo * Carlos Cordon-Cardo * Karla Paniagua Rivera * Mario Flores * Yu-Chiao Chiu * Yufei Huang * Shou-Jiang Gao ## Abstract Despite intensive studies during the last 3 years, the pathology and underlying molecular mechanism of coronavirus disease 2019 (COVID-19) remain poorly defined. Here, we examined postmortem COVID-19 lung tissues by spatial single-cell transcriptome analysis (SSCTA). We identified 18 major parenchymal and immune cell types, all of which are infected by SARS-CoV-2. Compared to the non-COVID-19 control, COVID-19 lungs have reduced alveolar cells (ACs), and increased innate and adaptive immune cells. Additionally, 19 differentially expressed genes in both infected and uninfected cells across the tissues mirror the altered cellular compositions. Spatial analysis of local infection rates revealed regions with high infection rates that are correlated with high cell densities (HIHD). The HIHD regions express high levels of SARS-CoV-2 entry-related factors including ACE2, FURIN, TMPRSS2, and NRP1, and co-localized with organizing pneumonia (OP) and lymphocytic and immune infiltration that have increased ACs and fibroblasts but decreased vascular endothelial cells and epithelial cells, echoing the tissue damage and wound healing processes. Sparse non- negative matrix factorization (SNMF) analysis of neighborhood cell type composition (NCTC) features identified 7 signatures that capture structure and immune niches in COVID-19 tissues. Trajectory inference based on immune niche signatures defined two pathological routes. Trajectory A progresses with primarily increased NK cells and granulocytes, likely reflecting the complication of microbial infections. Trajectory B is marked by increased HIHD and OP, possibly accounting for the increased immune infiltration. The OP regions are marked by high numbers of fibroblasts expressing extremely high levels of COL1A1 and COL1A2. Examination of single-cell RNA-seq data (scRNA-seq) from COVID-19 lung tissues and idiopathic pulmonary fibrosis (IPF) identified similar cell populations primarily consisting of myofibroblasts. Immunofluorescence staining revealed the activation of IL6-STAT3 and TGF-²-SMAD2/3 pathways in these cells, which likely mediate the upregulation of COL1A1 and COL1A2, and excessive fibrosis in the lung tissues. Together, this study provides an SSCTA atlas of cellular and molecular signatures of fatal COVID-19 lungs, revealing the complex spatial cellular heterogeneity, organization, and interactions that characterized the COVID-19 lung pathology. Key words * SSCTA * spatial single-cell transcriptome analysis * COVID-19 * coronavirus disease 2019 * SARS-CoV-2 * SNMF * sparse non-negative matrix factorization analysis * NCTC * neighborhood cell type composition analysis * trajectory inference * organizing pneumonia * fibrosis * IL6-STAT3 * TGF-β-SMAD2/3 ## Introduction Multiple single-cell RNA-seq (scRNA-seq) analyses of coronavirus disease 2019 (COVID-19) patients with different severities have improved our understanding of cellular diversity associated with infection and provided important molecular insights into the host immune response [1–5]. Despite intensive studies in the last 3 years, the pathology and underlying molecular mechanism of COVID-19 remain unclear. Severe COVID-19 is often accompanied by diffuse alveolar damage (DAD) that presents complex pathological manifestations and is heterogeneous within infected tissues and across patients [6]. The dissociation of tissue localization of cells in scRNA-seq has become a bottleneck to decoding the pathology features at the molecular and cellular levels and failed to reveal the immune signatures of the microenvironment for severe COVID-19. Spatial single-cell transcriptome analysis (SSCTA) promises to reveal the molecular basis of cellular heterogeneity, organization, and interactions in tissues and organs [7, 8]. However, analysis of these complex datasets including defining the spatial cellular organizations, immune microenvironment patterns, cell-cell interactions, and molecular signatures associated with disease pathophysiology remains a daunting task. Here, we utilized SSCTA to examine postmortem lung tissues from 5 cases with severe COVID-19 and one case without COVID-19. From 10,414,863 detected transcripts of cellular 221 genes from six tissues, we identified 1,719,459 cells that were mapped to 18 major parenchymal and immune cell types, all of which are infected by SARS-CoV-2. We further identified the spatial cellular and molecular signatures that define the patterns of SARS-CoV-2 infection, structural and pathological presentations, and associated immune microenvironments, which project the trajectories of disease progression. Together, this study provides an atlas of cellular and molecular signatures of fatal COVID-19 lungs and reveals the complex spatial cellular heterogeneity, organization, and interactions that characterize COVID-19 lung pathology. ## Results ### Spatial single-cell transcriptome analysis, cell segmentation, cell typing, and spatial mapping of cells Postmortem lung tissues from five COVID-19 autopsies and one postmortem case without COVID-19 were subjected to SSCTA (Fig. 1A). All five COVID-19 cases contracted SARS-CoV-2 in the first wave of the pandemic and had underlying conditions [9], including hypertension (cases 1-3), HIV infection and asthma (case 2), and Parkinson’s and chronic kidney diseases (cases 4 and 5). Hematoxylin and eosin (H&E) stains revealed various degrees of diffuse alveolar diseases (DAD), pulmonary thromboembolism, and lymphocytic infiltration in all cases (Fig. S1). Tissues 1 and 2 (1- 2C and 2-1A) had prominent organizing pneumonia (OP) or organizing diffuse alveolar damage while edema, hyaline membrane, and fibrin clot or microthrombi were prominent in tissues 3, 4, and 5 (3-1A, 4-3B and 5-3B). We designed probes for detecting the SARS-CoV-2 genome and 221 cellular genes covering markers of common lung parenchymal and immune cells, and immune and inflammatory genes induced by viral infections (Table S1). Following hybridization and gene decoding by *in- situ* sequencing, we stained the tissues with 4′,6-diamidino-2-phenylindole (DAPI) (Fig. S2) to facilitate cell segmentation (Fig. S2). ![Figure 1:](http://medrxiv.org/https://www.medrxiv.org/content/medrxiv/early/2023/02/26/2023.02.24.23286388/F1.medium.gif) [Figure 1:](http://medrxiv.org/content/early/2023/02/26/2023.02.24.23286388/F1) Figure 1: Spatial organization of parenchymal and immune cells in COVID-19 lung tissues revealed by spatial single-cell transcriptome analysis (SSCTA). A. Schematic illustration of the spatial single-cell sequencing pipeline used to generate the target dataset comprised of healthy control (non-COVID-19 tissue, n=1) and COVID-19 tissues (n=5). Samples were hybridized and decoded by sequencing. DAPI and H&E staining were carried out. All samples were annotated by an expert pathologist (Fig. S1). B. SSCTA workflow illustrating cell segmentation, cell typing, neighborhood cell- type composition (NCTC) analysis, and tissue pathology analysis. C. Average expression of cell type markers of all segmented cells across identified cell types in the SSCTA dataset. D. The UMAP projection of the gene expressions of segmented cells colored based on cell type, infection, and sample identifier. The UMAP shows negligible batch differences between the samples. E. Spatial visualization of identified cell types in two COVID-19 tissues and the non-COVID-19 tissue along with individual spatial plots illustrating the neighborhood cell type composition of specific cell types. We identified DAD in all COVID-19 tissues and orderly distribution of major parenchymal cells in non- COVID-19 tissue (PBC-PR). Further plots are provided in Fig. S1. F. Bar plots showing the percentages of identified cell types and their percentages of SARS-CoV-2 infected cells. A workflow was developed to analyze the SSCTA data (Fig. 1B). We detected 869,453 to 3,424,675 reads in the six lung tissues for a total of 10,414,863 reads (Table S2 and Fig. S3A). Using the Baysor algorithm [10], we segmented ∼89%-95% of reads into cells based on DAPI staining and the spatial distribution of reads (Fig. S4). We identified 186,659 to 470,294 cells for each sample for a total of 1,719,459 cells in these tissues (Table S2 and Fig. S3B). Over 99% of the cells harbored at least 5-15 reads (Fig. S5). We further filtered the cells deemed low quality (see Methods) and retained a total of 779,137 cells for the subsequent analyses. We summarized the reads for each gene and normalized the gene reads of segmented cells using *scTransform* [11]. We assigned cell types to 70-84% of the segmented cells based on the expressions of cell type-specific markers (Table S1) and identified a total of 18 cell types including 11 types of parenchymal cells and 7 types of immune cells in all six tissues (Fig. 1C and S6). We used Uniform Manifold Approximation and Projection (UMAP) to visualize the relationship of gene expressions of individual cells associated with different cell types in all tissues (Fig. 1D and S7). There were negligible batch differences as cells from different samples mixed well in individual cell clusters (Fig. 1D and S7). We noticed separated clusters for major parenchymal cells including alveolar cells (ACs) and fibroblasts. By contrast, immune cells were grouped together and mixed with vascular endothelial cells (VECs) in several clusters (Fig. 1D). The poor separation of immune cell types was largely due to highly expressed immunoglobulin kappa light chain (IGKC) and cathepsin L (CTSL) in these cell types (Fig. S8). We then mapped the individual cells to their spatial locations in the tissues (Fig. S9) and performed the neighborhood cell type composition (NCTC) analysis by computing a vector of 18 cell-type percentages in the neighborhood of individual cells. NCTC reveals spatial variations of the local cell type composition and informs local interplays among cell types. In agreement with the observed pathology that showed DAD in all COVID-19 tissues (Fig. S1), NCTC uncovered disorganized distributions of various cell types, especially the major parenchymal cells such as ACs, VECs, and fibroblasts (Fig. 1E, S1, and S10-S14). In contrast, organized structures and orderly distributions of these cell types were observed in the non-COVID-19 tissue (PBC-PR) (Fig. 1E, S1, and S15). Examination of regions with blood vessels revealed the lining of VECs along the vessels together with other vessel-associated cells including ACs, fibroblasts, and smooth muscle cells (SMCs) (Fig. 1E and S16). Furthermore, we observed abundant infiltrating immune cells including macrophage and monocytes (MMs), natural killer (NK) cells, and T- and B-cells along regions of blood vessels (Fig. 1E and S10-S14). We also observed the expected cell compositions of other structures including bronchiole, capillary, and endothelium (Fig. S17). These results validated the marker-based cell typing approach. ### SARS-CoV-2 infection alters the cell compositions of major parenchymal cells and induces immune infiltrations In agreement with the results of our previous study [9], we found SARS-CoV-2 reads in diverse cell types with infection rates ranging from 0.6% to 5% except that no infected cell was detected in the small number of erythrocytes (55) identified in tissue 1 (Fig. 1D, 1F, S7 and Table S3). Most infected cells were identified as ACs, fibroblasts, VECs, and MMs as they were the major cell types in the COVID-19 tissues (Fig. 1F). Consistent with the tissue pathology, we observed significant alterations in the compositions of different cell types in COVID-19 tissues compared to the non-COVID- 19 tissue (Fig. 1E and 1F). Among the parenchymal cells, there were reduced numbers of ACs from 44.08% to 12.81-28.23% and SMCs from 5.22% to 1.72-3.05%, and increased numbers of epithelial cells (ECs) from 2.08% to 2.40-8.42%, ionocytes from 0.22% to 0.41-1.36% and basal cells from 0.25% to 0.31-1.58%, suggesting their likely involvements in COVID-19 lung pathology (Fig. 1E, 1F, and Table S3). Furthermore, we observed increases in cell numbers of most types of the identified immune cells in COVID-19 tissues (Fig. 1E, 1F, and Table S3). In particular, there were increased inflammatory cells including MMs from 9.4% to 17.51-21.37% and NK cells from 1.00% to 2.32-9.19%, and adaptive immune cells including T cells from 2.40% to 4.34-6.81% and B cells from 1.91% to 2.70-10.51%, which were consistent with the reported inflammatory and cellular immune response following SARS-CoV-2 infection [12]. However, these changes were more subtle for some cell types in tissue 1, which had the least infected cells and mildest pathological manifestations with the best integrity of parenchymal cells among all COVID-19 tissues (Fig. 1E, 1F, S1, and S10). ### SARS-CoV-2 infection induces global differential gene expressions that marked pathological damages and inflammation with spatial cellular features We examined the differential gene expression between COVID-19 and non- COVID-19 lung tissues, and identified 9 upregulated and 10 downregulated genes, respectively (Fig. 2A, first panel). Comparison of SARS-CoV-2-infected cells in COVID- 19 tissues with all cells in the non-COVID-19 tissue confirmed the differential expressions in 16 of these 19 genes (Fig. 2A, second panel). Interestingly, differential gene expressions were also observed in 18 of these 19 genes between uninfected cells in COVID-19 tissues and all cells in the non-COVID-19 tissue (Fig. 2A, third panel). In contrast, there were only 3 differentially expressed genes between SARS-CoV-2- infected and -uninfected cells in the COVID-19 tissues (Fig. 2A, fourth panel). These results indicated that there were strong indirect effects such as those mediated by cytokines and complement activation induced by SARS-CoV-2 infection that likely contributed to the dysregulated gene expressions and pathology in COVID-19 lung tissues. Thus, we focused on comparing COVID-19 tissues with the non-COVID-19 tissue in subsequent analyses of dysregulated genes. ![Figure 2:](http://medrxiv.org/https://www.medrxiv.org/content/medrxiv/early/2023/02/26/2023.02.24.23286388/F2.medium.gif) [Figure 2:](http://medrxiv.org/content/early/2023/02/26/2023.02.24.23286388/F2) Figure 2: SARS-CoV-2 infection induces global differential gene expression that mediates pathological damages and inflammation with cell type and spatial features. A. Differential gene expressions between the identified cells in the COVID-19 non-COVID-19 tissue with various conditions are illustrated. The title indicates the conducted differential expression analysis. Within each row, the bubble size indicates the -log10 of the corrected p-value and the color indicates the log2 fold change of the corresponding gene. Each column indicates the different COVID-19 tissues. The first three panels indicate differential expression analysis between COVID-19 tissues and non-COVID-19 tissue while the fourth panel indicates the differential expression analysis between the infected and uninfected cells in the respective COVID-19 tissues. B. Spatial gene expression map of highly expressed CTSL and IGKC. C. Spatial gene expression maps of the rest of the differentially expressed genes. D. Pie chart illustrating the percentages of cell types in which the genes are expressed in each tissue. E. Differential expression analysis of COVID-19 tissues compared to non- COVID-19 tissue in cell types that had significant changes in cell numbers based on Fig. 1F. The expressions of the 19 dysregulated genes in different cell types and their tissue distributions had significant variations suggesting their complex involvements in different aspects of SARS-CoV-2 infection and COVID-19 lung pathology (Fig. 2B, 2C, S8, and S18-S23). We determined the contributions of individual cell types to the expressions of these genes across all the tissues (Fig. 2D). Of the 19 genes, 10 of them were cell markers (Table S1). All 6 parenchymal cell markers including SFTPC, SFTPA1, and advanced glycosylation end product (AGE) receptor (AGER) for ACs, INMT for fibroblasts, smooth muscle actin alpha 2 (ACTA2) for SMCs, and claudin 5 (CLDN5) for VECs were downregulated, and mainly expressed in their respective cell types (Fig. 2A, first panel, 2C, 2D, S8 and S18-S23), of which ACs and SMCs had decreased cell numbers in COVID-19 tissues (Fig. 1E, 1F and S7). Both SFTPC and SFTPA1 are implicated in lung homeostasis and functions, and their mutations and dysregulation are associated with pulmonary fibrosis [13–15]. Additionally, SFTPA1, a C-type lectin, which binds to specific carbohydrate moieties on lipids and the surface of microorganisms, is essential in the defense against respiratory pathogens [16–18] while the pulmonary-associated surfactant protein C encoded by SFTPC is a marker for COVID-19 patients with high viral loads [19]. AGER, predominantly expressed in ACs (Fig. 2C), is a multiligand receptor interacting with AGE and other molecules implicated in lung homeostasis, development and inflammation, and certain diseases such as diabetes and Alzheimer’s disease, and regulates diverse pathways including MyD88-dependent, nuclear receptors, TNF-α, ERK1/2 and p38 MAPK, and p53/TP53 pathways [20–23]. Following interaction with S100A12, AGER triggers the activation of mononuclear phagocytes, lymphocytes, and endothelium by generating key pro- inflammatory mediators [24]. Indeed, AGER-related pathways are activated by SARS- CoV-2 infection and are implicated in the COVID-19 lung pathology [25, 26]. ACTA2 encodes a smooth-muscle actin involved in lung functions including cell motility, tissue structure and integrity, and intercellular signaling [27]. Dysregulation of ACTA2 is linked to a variety of vascular diseases as well as multisystemic smooth muscle dysfunction syndrome [28, 29]. CLDN5 is an integral membrane protein and component of tight junction strands regulating the integrity of epithelial or endothelial cell sheets, and immune cell transmigration [30–32]. SARS-CoV-2 infection suppressed CLDN5 expression contributing to the disruption of respiratory vascular barriers [33]. Thus, downregulations of parenchymal cell markers including SFTPC, SFTPA1 and AGER in ACs, ACTA2 in SMCs, and CLDN5 in VECs might be involved in SARS-CoV-2 infection and COVID-19 lung pathology (Fig. 2D). Several other downregulated genes including caveolin 1 (CAV1) and alpha-2- macroglobulin (A2M) are implicated in numerous lung functions (Fig. 2A). CAV1, a component of caveolae involved in multiple pathways such as integrins-mediated and Ras-ERK signaling, is essential for lung functions and during the host response to infections [34, 35]. A2M is an inhibitor of broad spectrum of proteases including trypsin, thrombin, and collagenase [36, 37]. It disrupts inflammatory cascades by inhibiting inflammatory cytokines, is implicated in Alzheimer’s disease, and regulates extracellular matrix organization and platelet cytosolic Ca2+ [36, 37]. CAV1 and A2M were downregulated and mainly expressed by ACs, fibroblasts, VECs, and MMs (Fig. 2A, first panel, 2D, S8, and S18-S23), of which ACs and fibroblasts had decreased while MMs had increased cell numbers in COVID-19 tissues (Fig. 1F). Three upregulated immune marker genes including complement component C1q, A chain (C1QA) and B chain (C1QB), and neural cell adhesion molecule 1 (NCAM1) are involved in innate immune response and inflammation (Fig. 2A, first panel) [38, 39]. C1QA and C1QB were predominantly expressed in MMs (Fig. 2D, S8 and S18-S23). SARS-CoV-2 infection induced robust complement activation, contributing to tissue damage and COVID-19 lung pathology [40–42]. Both cathepsin L (CTSL) and granzyme K (GZMK) could be involved in complement activation [43, 44]. CTSL was induced after SARS-CoV-2 infection and enhanced SARS-CoV-2 infection [45] while GZMK, a serine protease from the cytoplasmic granules of cytotoxic lymphocytes (CTL) and NK cells that recognize and lyse specific target cells [46], could limit the spread of SARS-Co-V-2. Both CTSL and GZMK were upregulated and mainly expressed by ACs, fibroblasts, MMs, and VECs while CTSL exhibited a comparable expression level but with extremely high levels in localized regions in COVID-19 tissues (Fig. 2A, 2B, 2C, 2D, S8 and S18-S23). NCAM1, the marker of NK cells, is a member of the immunoglobulin superfamily involved in cell-to-cell interactions as well as cell-matrix interactions [47]. The upregulations of C1QA, C1QB, CTSL, GZMK, and NCAM1 were consistent with the increased numbers of MMs and NK cells in COVID-19 tissues (Fig. 1F, S8, and S18- S23). Thus, increased infiltrations of MMs and NK cells were likely involved in COVID- 19 lung complement activation and inflammation. Two immune-modulating cytokines were upregulated including colony-stimulating factor 3 (CSF3), a member of the IL-6 superfamily of cytokines that controls the production, differentiation, and function of granulocytes involved in the innate immune response [48, 49], and inducible T cell costimulator ligand (ICOSLG or B7-H2) involved in positive regulation of interleukin-4 production and CD28 signaling in T-helper cells [50, 51]. The major cells that expressed CSF3 and its receptor CSF3R, as well as ICOSLG and its receptor ICOS, were ACs, fibroblasts, VECs, and MMs (Fig. 2D). Numerous other upregulated genes in COVID-19 tissues are involved in inflammatory and immune responses (Fig. 2A). Leukocyte-specific transcript 1 (LST1), a myeloid leukocyte-specific transmembrane adaptor protein recruiting protein tyrosine phosphatases SHP-1 and SHP-2 to the plasma membrane, inhibits the proliferation of lymphocytes, and its expression is enhanced by lipopolysaccharide (LPS), interferon- gamma (IFN-ψ) and bacteria [52, 53]. Immunoglobulin kappa light chain (IGKC) is essential for antibody production but is also expressed in non-lymphoid cells [54–56]. Both LST1 and IGKC were mainly expressed in ACs, fibroblasts, VECs and MMs (Fig. 2B-E, S8 and S18-S23), which might mediate dysregulation of immune cells in COVID- 19 lung tissues. Finally, fibulin 1 (FBLN1), a secreted glycoprotein incorporated into a fibrillar extracellular matrix in a calcium-dependent manner and mediates platelet adhesion via binding fibrinogen [57, 58], was mainly expressed in ACs, fibroblasts, VECs, and MMs (Fig. 2D). We examined differential gene expressions in cell types that had significant changes in cell numbers (Fig. 1F). Compared to non-COVID-19 tissue, ACs from COVID-19 tissues had 20 differentially expressed genes, of which 14 were identified at the whole tissue level (Fig. 2A, first panel and 2E). The most upregulated genes were LST1, CTSL, IGKC, ICOSLG, C1QB, GZMK, and NCAM1 while the most downregulated genes were SFTPC, A2M, SFTPA1, CLDN5, and CAV1 (Fig. 2E), which likely contribute to the decreased cell numbers in COVID-19 tissues (Fig. 1F). Several differentially expressed genes were also identified in fibroblasts, VECs and ECs (Fig. 2E). Among the immune cells, MMs had the most differentially expressed genes including upregulation of C1QB, NCAM1, and IGKC, and downregulation of A2M, CD68, TYROBP and PRG4 (Fig. 2E). ### Spatial analysis identifies regions with high SARS-CoV-2 infection rates that match high cell densities, high levels of viral entry-related factors, and localized pathology Even though we did not observe notable global differences in gene expression between infected and uninfected cells within each COVID-19 tissue, the spatial nature of the SSCTA permitted the identification of the local impact of SARS-CoV-2 infection on lung pathology. We performed the local infection rate analysis to examine the spatial distributions of local SARS-CoV-2 infection (Fig. 3A) and identified regions with high and low infection rates in each tissue (Fig 3B; local Moran’s I, p-value < 0.05). SARS- CoV-2 infects host cells by binding to its entry receptor, angiotensin-converting enzyme 2 (ACE2), and subsequent engagement of host proteases and other entry-related factors including FURIN, transmembrane protease, serine 2 (TMPRSS2), and neuropilin 1 (NRP1) [59]. We assessed the spatial correlation between high-infection regions and cells expressing transcripts of these viral entry-related proteins. We found a significant correlation between high infection rates and high expressions of ACE2, FURIN, TMPRSS2, and NRP1 (Fig. 3C; bivariate local Moran’s I, p-value < 0.05), supporting the essential roles of these cellular proteins in SARS-CoV-2 infection in lung tissues. ![Figure 3:](http://medrxiv.org/https://www.medrxiv.org/content/medrxiv/early/2023/02/26/2023.02.24.23286388/F3.medium.gif) [Figure 3:](http://medrxiv.org/content/early/2023/02/26/2023.02.24.23286388/F3) Figure 3: Spatial analysis of local SARS-CoV-2 infection rates identifies regions with high infection rates that match high cell densities, high levels of viral entry- related factors, and localized pathology. A. Spatial distribution of the SARS-CoV-2 infected cells. The spatial visualization of the infection rates revealed infection hotspots in the COVID-19 tissues. B. Moran’s I analysis of spatial patterns revealed spatially distinct high and low infection regions in each COVID-19 tissue. C. Spatial visualization of the significant local regions identified by Bivariate Moran’s I analysis, in which infection rates and the expressions of TMPRSS2, NRP1, FURIN and ACE2 are spatially correlated. D. Spatial visualization of the kernel density plot for the cell densities revealed highly dense cellular hotspots in the COVID-19 samples. E. Spatial visualization of the significant local regions identified by Bivariate Moran’s I analysis, in which cell densities and the infection regions are spatially correlated, identified as high- density high-infection (HIHD) regions. F. Spatial visualization of the ligand-receptor coexpression maps for the CSF3-CSF3R pair, which are highly expressed in the high infection regions. G. Spatial visualization of the significant local regions identified by Bivariate Moran’s I analysis, in which both high infection and high CSF3-CSF3R coexpression regions are spatially correlated. Previous studies reported increased densities of both parenchymal and immune cells in COVID-19 lung tissues [6]. Here, we observed a positive spatial correlation between cell densities and SARS-CoV-2 infection rates especially in tissues 1-4 (Fig. 3D, Table S4), suggesting that increased cell densities might be associated with high infection rates. We segmented regions with high infection rates and high cell densities (HIHD) (Fig. 3E; bivariate local Moran’s I, p-value < 0.05). Compared with H&E staining (Fig. S1), we found a close association between HIHD and OP, which was most prominent in tissues 1 and 2, as well as lymphocytic and immune infiltration, which was most severe in tissues 3 and 4 (Fig. 3E). Both are the most common lung pathological manifestations in COVID-19 lungs [12]. Hence, high infection rates might cause persistent damages to the lung tissues, resulting in the infiltration of immune cells, and increased wound repair and fibrosis [60, 61]. We also observed HIHD areas adjacent to blood vessels in tissues 1 and 3 (Fig. 3E and S1), which might be due to the homing of the circulating virus and immune cells in the bloodstream to these sites or that these cells were more susceptible to SARS-CoV-2 infection [9, 62]. ### Regions with high SARS-CoV-2 infection rates have unique cellular and gene expression features We next examined cell type compositions of the high and low SARS-CoV-2 infection regions. Compared to low-infection regions across the tissues, the compositions in high-infection regions exhibited significant increases in fibroblasts, echoing the observed OP in HIHD areas, and an increase in ACs, possibly due to the proliferation of ACs response to injury repair of alveoli and alveoli capillaries (Fig. 3E and 4A). However, there were decreases in ECs and VECs, which might reflect the damages caused by the infection (Fig. 3E and 4A). These results suggested that high- infection regions might suffer more damages. However, except for T cells and granulocytes, the compositions of major immune cell types including MMs between high and low infection regions remained similar (Fig. 4A), suggesting induction of a broad immune infiltration by SARS-CoV-2 infection, as well as the observed broad impact on the tissue-wide gene expression (Fig. 2A). ![Figure 4:](http://medrxiv.org/https://www.medrxiv.org/content/medrxiv/early/2023/02/26/2023.02.24.23286388/F4.medium.gif) [Figure 4:](http://medrxiv.org/content/early/2023/02/26/2023.02.24.23286388/F4) Figure 4: Spatial analysis reveals cellular and gene expression features in regions with high SARS-CoV-2 infection rates. A. A bar plot showing the cell type composition of all identified cell types in the high and low infection regions. The increase in alveolar cells (ACs) might suggest the proliferation of ACs due to injury repair of alveoli and alveoli capillaries and a decrease in epithelial cells (ECs) and vascular endothelial cells (VECs) might suggest the damages caused by the infection. **B.** UMAP projection of the NCTC analysis results of high infection regions (refer to figure 3B) revealed three distinct clusters out of which two were distinctly enriched in local compositions of ACs and fibroblasts, which are illustrated on the right with their spatial visualizations. **C.** Spatial visualization of the compositions in high infection regions revealed negative correlation of fibroblasts and ACs (-0.6, Pearson correlation) which overlapped with annotated OP regions (high fibroblasts, low ACs). **D.** Differential expressions of high and low infection regions of the COVID-19 tissues compared to the non-COVID-19 tissue, which share 18 of the 20 differentially expressed genes, confirming strong indirect effect of SARS-CoV-2 infection. **E.** Differential expressions between the high and low infection regions of the COVID-19 tissues revealed distinct differentially expressed genes showing upregulation of AC markers SFTPA1 and SFTPA2 in tissues 1, 3, and 5, and fibroblast markers COL1A1 or COL1A2 in all COVID-19 tissues. **F.** Spatial visualization of the spatial gene expression map of COL1A1 and COL1A2 markers. The patchwork of high infection regions in tissue 1 prompted us to inspect patterns of local cell-type compositions. Interestingly, NCTC analysis showed three clusters in high-infection regions, of which two were distinctly enriched with prominent local compositions of fibroblasts and ACs, respectively (Fig. 4B). Examining their spatial distributions confirmed a negative *in-situ* correlation between local compositions of fibroblasts and ACs (-0.6, Pearson correlation), with regions showing high (low) compositions of fibroblasts accompanied with low (high) compositions of ACs (Fig. 4B and C). Strikingly, the local regions with high fibroblasts but low ACs were mostly annotated OP regions (Fig. S1). We next investigated the gene expression patterns in the high- and low-infection regions. Compared with the normal tissue, the high and low infection regions presented highly consistent differential expression patterns as they shared 18 of the 20 differentially expressed genes (Fig. 4D), confirming the strong indirect effect of SARS- CoV-2 infection (Fig. 2A). One of the two uniquely differentially expressed genes in high infection regions was again cytokine CSF3, which was upregulated in all infected tissues. We, therefore, investigated spatial co-expression of CSF3 and its receptor CSF3R (Fig. 3F) and found that their co-expression patterns also presented spatial correlations in high infection regions (Fig. 3G; bivariate local Moran’s I, p-value < 0.05), suggesting that the CSF3-CSF3R axis might modulate the immune response in high- infection regions. In contrast to the consistent differentially expressed patterns of both high- and low-infection regions when compared with the normal control, direct comparison between high- and low-infection regions revealed distinct differentially expressed patterns. Particularly, we observed, in high- *vs.* low-infection regions, upregulation of AC markers SFTPA1 or SFTPA2 in tissues 1, 3, and 5 (Fig. 4E and 2C) and fibroblast markers COL1A1 or COL1A2 in all COVID-19 tissues (Fig. 4E and 4F), despite that SFTPA1 and COL1A2 displayed down-regulation in either high- or low- infection regions in COVID-19 tissues when compared to the normal tissue (Fig. 4D). The increased expressions of these AC and fibroblast markers were also consistent with the observed increases in ACs and fibroblasts in high-infection regions in the corresponding tissues (Fig. 4A). Cell type-specific differential expressions further confirmed the upregulation of SFTPA1 or SFTPA2 in ACs, and COL1A2 in fibroblasts in infected tissues (Fig. 4E and S24). However, different cell types exhibited complex differential expression patterns between high *vs.* low infection regions in individual COVID-19 tissues (Fig. S24), likely due to the distinct COVID-19 pathologies or stages associated with the individual infected tissues. ### Sparse non-negative matrix factorization analysis identifies seven cell composition signatures that recapitulate different healthy and disease statuses Although different DAD phases can inform COVID-19 severity, they are confounded by many factors including the underlying conditions of the patients as illustrated by our five COVID-19 cases (Fig. S1). Importantly, the complex tissue microenvironments of parenchymal and immune cells associated with different stages of DAD within and across tissues have not been well characterized. Therefore, we determined the extent to which tissue-independent and tissue-specific spatial patterns in spatially resolved cell-type compositions are associated with different DAD stages. For this purpose, we performed the sparse non-negative matrix factorization (SNMF) [63] of the NCTC vectors of all cells in all tissues. SNMF decomposed these NCTC vectors into a sparse linear combination of cell-type composition signatures that define spatial cell- type patterns in these tissues. We obtained a factorization of seven NCTC signatures after assessing the compactness and biological meanings of the factorization results for different maximum numbers of factors (Fig. 5). Close examination of cell type compositions of the signatures (Fig. 5A and Table 5) and their spatial prevalence across all tissues (Fig. 5B) defined by the factor loading revealed their associations with normal structures, and broad and tissue-specific infection. Specifically, Signature 1 likely defined the cell type composition of “normal-like lung alveoli” because it has 56.4% ACs, 18.5% VECs, 11.9% fibroblasts, and 5.3% SMs, which were comparable to the percentages of the non-COVID-19 tissue (Fig. 1F and Table S3). Thus, the high-loading regions of Signature 1 in COVID-19 tissues might mark the less damaged regions. ![Figure 5:](http://medrxiv.org/https://www.medrxiv.org/content/medrxiv/early/2023/02/26/2023.02.24.23286388/F5.medium.gif) [Figure 5:](http://medrxiv.org/content/early/2023/02/26/2023.02.24.23286388/F5) Figure 5: Sparse non-negative matrix factorization (SNMF) identifies seven cell composition signatures that recapitulate different healthy and disease statuses. **A.** SNMF revealed seven NCTC signatures which are illustrated in a logo plot. **B.** Spatial NCTC signatures visualized after factor loading for each tissue revealed their associations with normal structures and broad and tissue-specific infection. **C.** Trajectory analysis of the cell type compositions with Signatures 1 and 3 generated a trajectory of reduced normal alveoli and increased immune infiltration from the non- COVID-19 tissue to COVID-19 tissues 1, 4, 3, 5, and 2. The trajectory agrees with the order of the severity of pathology in the tissues identified in the H&E stained images (Fig. S1). Naturally, the non-COVID-19 tissue had a high prevalence of Signature 1 (Fig. 5B and S25). Tissue 1 had the highest prevalence of Signature 1 among all COVID-19 tissues reflecting its least damaged pathology, followed by tissue 4, 3, 2, and 5 (Fig. 5B, S1, and S25). Signature 2 described a composition of 61% ECs, 14% ciliated cells, and 10% basal cells. It had sparse but high-intensity loadings in localized regions identified as bronchial tubes (Fig. 5B and S1). In contrast to Signatures 1 and 2 which characterized cell-type compositions of normal structures, Signature 3 captured the “broad immune infiltration” as it contained 28.6% B cells, 22.3% MMs, 17.5% T cells, 9.6% granulocytes, and 8.7% ECs (Fig. 5A and Table S5). As expected, it was the least common in the non-COVID-19 tissue (Fig. 5B and S25). Its loading distributions in COVID-19 tissues were consistent with the lowest immune infiltration in tissue 1, followed by 4, 3, 2, and 5 (Fig. 5B and S25), closely corroborating the order of damage revealed by Signature 1. Trajectory analysis of the cell type compositions with Signatures 1 and 3 using *Slingshot* [64] defined a trajectory of reduced normal alveoli and increased immune infiltration from the non-COVID-19 tissue to tissue 1, 4, 3, 5, and 2 (Fig. 5C), which was in agreement with the order of their increased severity of pathology observed in H&E staining (Fig. S1). Remarkably, Signatures 4 to 7 captured four distinct types of immune microenvironments or niches in these COVID-19 tissues. Signature 4 was characterized by high ACs (64.8%) with MMs (27.3%), likely reflecting the infiltration of alveolar MMs and proliferation of ACs in wound healing, all of which are features of exudative DAD [65]. The relatively higher prevalence of Signature 4 in the less severe tissues 1 and 4 suggested its association with early phases of DAD (Fig. 5B and S25). There were also overlaps in regions with higher loadings of Signature 4 with the annotated regions of immune infiltration in the H&E image in tissue 4 (Fig. 5B and S1). Signatures 5 and 6 depicted inflammations around VECs, where Signature 5 was associated with extensive infiltration of MMs (46.7%) and T cells (8.7%) around VEC (26.3%), whereas Signature 6 was associated with substantial infiltration of NK cells (31.8%) and granulocytes (13.6%) around VECs (27.6%) and ECs (12.2%). These results were consistent with the influx of MMs and NK cells (Fig. 1E and F), which was likely due to the activation of complement and coagulation systems (Fig. 2A and C), resulting in damages and therefore reduced numbers of ACs and VECs (Fig. 1E and F). Signature 5 was more common in the less severe tissues 1, 3 and 4 but Signature 6 appeared predominantly in 3 (Fig. 5B and S25), particularly in low infection region (Fig. 3A). Indeed, lymphocytic infiltrations were also identified in the regions associated with Signatures 5 and 6 (Fig. S1). Signature 7 denotes extremely high fibroblasts (73.4%) with infiltration of MMs (15.8%). Advanced disease in COVID-19 patients often has an excessive accumulation of fibroblasts and may develop pathological fibrosis due to chronic inflammation from the infiltrated immune cells and dysregulation of the extracellular matrix triggered (ECM) remodeling [66]. As expected, there was a high correlation between the COL1A1 expression level (Fig. 4F) with Signature 7. The loadings of this signature were high in tissues 1 and 2, which were annotated as OP regions, and in tissue 5, which were annotated as fibrin and hyaline membrane in addition to OP regions (Fig. 5B, S1 and S25), suggesting a connection with OP, fibrin and hyaline membrane. ### Unique cell composition signatures define two separate pathological trajectories We then applied Signatures 4 to 7 to delineate the trajectories of DAD characterized by the immune responses against viral infection. We added Signature 1 “normal-like alveoli” so that the non-COVID-19 tissue could be included as a reference. We used these 5 signatures to redefine the landscape of spatially resolved cell type compositions in these tissues. The UMAP visualization (Fig. 6A) recapitulated the tissue severity inferred by Signatures 1 and 3 (Fig. 5C) but presented more complex patterns, with each immune signature defining a largely distinct group of cell populations (Fig. 6B). We applied *Slingshot* to infer pseudo-progressions of cells defined by these signatures and identified two trajectories (Fig. 6A). Clustering of cells based on pseudo- time values further segmented the trajectories into 6 stages (Fig. 6A). At the onset, the two trajectories shared a common path including stages T1 and most T2, which traversed the cells enriched by Signature 1 “normal-like alveoli” and Signature 4 “ACs with MMs infiltration” in T1 and then those enriched by Signature 4 and Signature 5 “VECs with MMs infiltration” in T2 (Fig. 6A). These stages were consistent with the less damaged presentations of tissues 1 and 4 (Fig. 6D and S1) and showed a reduction of ACs due to likely lung injury and an increase of MMs possibly promoted by innate immune responses (Fig. 6C). Inspecting associated tissue regions in the H&E images showed characteristics of early DAD manifested by shedding of ACs, capillary leakage, lymphocytic/immune infiltration, and yet preserved alveoli structure (Fig. 6D and E, T1 and T2). After the shared section, the trajectories diverged into two paths. The first path (Trajectory A) included T3a and T4a and progressed through cells mostly from tissues 3 and 4 enriched in Signature 5 “VEC with MMs infiltration” ending in Signature 6 “VEC with NK cell infiltration”, which was almost exclusively from tissue 3 (Fig. 6A, B and D). We observed a persistent decrease in ACs but an increase in VECs, NK cells and B cells (Fig. 6C). The H&E staining of tissue regions for T3a and T4a showed hyaline membrane, fibrin, lymphocytic infiltration, and collapsed alveoli and capillary structure (Fig. 6E, T3a and T4a, and S1). These features are consistent with a more advanced DAD. The close association of this path with tissues 3 and 4 suggested that they encompassed similar DAD-related tissue patterns and immune niches. Indeed, the postmortem report showed similar findings for these two patients including bilateral pulmonary consolidation [9]. ![Figure 6:](http://medrxiv.org/https://www.medrxiv.org/content/medrxiv/early/2023/02/26/2023.02.24.23286388/F6.medium.gif) [Figure 6:](http://medrxiv.org/content/early/2023/02/26/2023.02.24.23286388/F6) Figure 6: Unique cell composition signatures define two separate pathological trajectories. A. The UMAP projection of the Signatures 1 (for non-COVID-19 tissue reference) and Signatures 4-7 recapitulated the severity of pathological progression identified in figure 5C. Two trajectories were derived from the pseudo-progressions of cells defined by these signatures. The trajectories showed a common shared onset T1 and T2 which traverses through Signature 1 “normal-like-alveoli” and Signature 4 “ACs with MMs infiltration” in T1 and then those enriched by Signature 4 and Signature 5 “VECs with MMs infiltration” in T2. The trajectory then diverged into two paths T3a and T4a, and T3b and T4b. Path/Trajectory A defined by T3a and T4a progressed through cells mostly from tissues 3 and 4 enriched in Signature 5 “VEC with MMs infiltration” ending in Signature 6 “VEC with NK cell infiltration”, which was almost exclusively from tissue 3. Path/Trajectory B defined by T3b and T4b are enriched by Signature 6 in T3b and Signature 7 “high fibroblasts with MMs infiltration” in T4b. **B.** UMAP projections of the five signatures colored based on their individual signature loadings. **C.** Line plots illustrating the cell type percentages in band T1, T2, T3(a, b), and T4(a, b). **D.** Spatial visualization of the cells by cell type in each band is shown along with an identified representative ROI (red box with dotted lines) for each COVID-19 tissue. Column 7 illustrates the H&E morphology image of the COVID-19 tissues with black dotted lines highlighting the high infection regions. **E.** Zoomed-in images of the ROIs selected are shown. A close examination revealed characteristics of early DAD in T1 and T2, hyaline membrane, fibrin, lymphocytic infiltration, collapsed alveoli and capillary structure in T3a and T4a, and edema, increased fibrin, and organizing pneumonia (OP) in T3b and T4b suggesting a more advanced stage of DAD. The second path (Trajectory B) included two stages, i.e., T3b and T4b, which are enriched by Signature 6 in T3b and Signature 7 “high fibroblasts with MMs infiltration” in T4b highlighted by increased fibroblasts and immune infiltration of especially MMs (Fig. 6A, B, and C). Their associated tissue regions were mostly in tissues 1, 2, and 5 (Fig. 6D), suggesting potentially similar DAD tissue patterns and underlying immune niches in these tissues. The postmortem findings indeed revealed massive pulmonary emboli that were the main cause of death for cases 1 and 2 [9]. Close examination of these regions in the H&E staining showed edema, increased fibrin, and OP (Fig. 6D, E, T3b and T4b, and S1), which are pathological signatures of DAD. There was a striking association between regions in tissues 1, 2 and 5 in T4b with the identified OP, fibrin and hyaline membrane regions and a significant positive correlation between the progression of Trajectory B and cell density (0.504, Spearman correlation, *p-value* < 0.05) suggesting a more advanced stage of T4b than T3b with tissue 1 manifested in a more localized while tissues 2 and 5 in a more systematic manner. These results demonstrated a high degree of inter and intra-tissue heterogeneity in DAD-associated patterns. While tissue 1 contained regions affiliated with both trajectories, tissues 2 and 5 were mostly associated with the trajectory featured by high fibroblasts and their reorganization, and tissues 3 and 4 were with the trajectory featured by increased VEC and immune infiltration. Next, we investigated the genes whose expressions in a cell type were correlated with the trajectories. We found 53 and 67 genes significantly correlated with Trajectory A and B in at least one cell type, respectively (Spearman correlation, *p-value*<0.05, permutation test; Table S6). In both trajectories, we obtained a negative correlation for ACs marker SFTPC in ACs (Fig. S26A, S26B, Table S6, and Table S7), suggesting a decrease in its expression, which echoes the reductions in ACs in both trajectories (Fig. 6C). We also found a positive correlation of NCAM1 in NK cells with Trajectory A and positive correlations of COL1A1 in Fibroblasts and C1QB in MMs with Trajectory B (Fig. S26A, S26B, Table S6, and Table S7), all of which were consistent with the observed changes of cell type compositions and DAD pathologies along the trajectories (Fig. 6C). To gain further insight into potential ligand-receptor interactions that underlie the immune patterns associated with the trajectories, we correlated the ligand-receptor co- expressions with the trajectories. We found that CD40LG-CD40, CD80-CD28 and CXCL10-CXCR3 had the highest positive correlations with Trajectory A in all key cell types including especially VECs, ACs and Fibroblasts (Fig. S26C) while CD40LG-CD40, CSF3-CSF3R and CXCL10-CXCR3 were the top positive correlated pairs with Trajectory B in ACs, Fibroblasts, VECs and MMs (Fig. S26D). Thus, these ligand- receptor pairs might contribute to the chemoattraction of the immune cells as a result of SARS-CoV-2 infection. ### IL6-STAT3 and TGF-β-SMAD2/3 pathways mediate COVID-19 lung fibrosis and organizing pneumonia Because Trajectory B progressed into tissue regions that manifested extensive OP, we investigated the spatial patterns and molecular signatures associated with OP, which was mostly observed in tissues 1 and 2 and to a less extent in tissue 5, which also contained fibrin and hyaline membrane (Fig. 7A and S1). Analysis of cell compositions in regions with OP revealed an overall increase in cell density, and obvious decreases in cell numbers of ACs and VECs but an increase in cell numbers of fibroblasts compared to other regions in the same tissues and the non-COVID-19 tissue (Fig. 7B, C and D). As in the whole tissues, we observed downregulations of AGER, CAV1, SFTPA1, SFTPC, and A2M, and upregulation of IGKC in OP regions compared to the normal tissue (Fig. 7E). The changes of AGER, SFTPA1, SFTPC, and A2M were primarily observed in fibroblasts, MMs and VECs with downregulation of A2M also being observed in ACs (Fig. 7F, G and H). The most striking observation is the upregulation of genes encoding for collagen type I alpha 1 and 2 chains (COL1A1 and COL1A2) in OP in both tissues compared to the normal tissue or to regions without OP in the same COVID-19 lung tissues (Figure 7E and G). Furthermore, both COL1A1 and COL1A2 were upregulated in fibroblasts in OP regions compared to fibroblasts from either normal tissue or non-OP regions from the same tissue (Fig. 7G, H and I). COL1A1 and COL1A2 form the triple helix of a fibril-forming collagen, which is essential for lung homeostasis. Dysregulation of COL1A1 and COL1A2 leads to lung inflammation and fibrosis [67, 68]. Both COL1A1 and COL1A2 were mainly expressed in fibroblasts (Fig. 7G and 7H), whose cell numbers were doubled in OP regions in both tissues compared to non-OP regions (Fig. 7D). Interestingly, the expression levels of COL1A1 and COL1A2 were highly heterogeneous in the individual cells in OP regions with extremely higher levels in 10-20% of the fibroblasts, which were absent in fibroblasts in non-OP regions or any other types of cells in the same COVID-19 tissues (Fig. 7J). COL1A1 and COL1A2 were not expressed in high levels in any types of cells in normal lung tissue (Fig. 7J). These results indicated that COVID-19 OP regions had severe dysregulation of fibroblasts with upregulation of COL1A1 and COL1A2, and abnormal ACs and VECs as a result of downregulation of AGER, SFTPA1, and SFTPC. ![Figure 7:](http://medrxiv.org/https://www.medrxiv.org/content/medrxiv/early/2023/02/26/2023.02.24.23286388/F7.medium.gif) [Figure 7:](http://medrxiv.org/content/early/2023/02/26/2023.02.24.23286388/F7) Figure 7: IL6-STAT3 and TGF-β-SMAD2/3 pathways mediate COVID-19 lung fibrosis and organizing pneumonia. A. Zoomed in OP region (bottom, red color border) and other non-OP regions (top, black color border) from the H&E image of COVID-19 tissue 1-2C are illustrated. **B.** Identified cells colored based on the cell types revealed differences in cell composition signatures between OP and other regions. **C.** Bar plot showing an increase in cell density in OP regions compared to other regions. **D.** Bar plot showing the percentage of cells per cell type in COVID-19 and non-COVID-19 tissues. There is a decrease in ACs and VECs, and an increase in fibroblasts in the OP regions of the COVID-19 tissue. An obvious increase in ACs is seen in the non-COVID- 19 sample. **E.** Differential gene expressions between the identified cells in the OP regions of the COVID-19 tissue sample with cells from the non-COVID-19 tissue and cells from the other non-OP regions of the same COVID-19 tissue is illustrated. We observed an upregulation of genes encoding for collagen type I alpha 1 and 2 chains (COL1A1 and COL1A2) in OP in both tissues compared to the normal tissue or to regions without OP in the same COVID-19 lung tissues. **F.** Pie chart illustrating the percentages of cell types in which the selected genes are expressed in tissues 1 and 2. **G.** Spatial visualization of differentially expressed genes from figure 7E. The size of the dots indicates the expression level, and the color indicates the cell type of the cell where the gene is expressed. **H.** Differential gene expressions between the identified cells in the OP regions of the COVID-19 tissue with cells from the non-COVID-19 tissue as well as non-OP regions from the same COVID-19 tissue based on their cell types. Both COL1A1 and COL1A2 were upregulated in fibroblasts in OP regions compared to fibroblasts from either normal tissue or non-OP regions from the same tissue. **I.** The density plot of the COL1A1 and COL1A2 expressions in OP and other non-OP regions in COVID-19 tissue samples 1-2C, 2-1A, and non-COVID-19 tissue is illustrated. **J.** Examination of the networks associated with the Pulmonary Fibrosis Idiopathic Signaling Pathways reported by the Ingenuity Pathway Analysis is illustrated. Both IL6-STAT3 and TGF-β-SMAD2/3 pathways mediate the expression of COL1A1 and COL1A2. **K.** Multi-color IFA staining revealed that cells with a high COL1A1 expression level also had high levels of IL6Rα, p-STAT3, TGF-βR2 and SMAD2/3, confirming the activation of these pathways and their potential roles in the upregulation of COL1A1. To further confirm that increased expressions of COL1A1 and COL1A2 were hallmarks of COVID-19 OP, we examined their expressions in different cell types in single-cell RNA-seq (scRNA-seq) data from COVID-19 lung tissues [4]. Consistent with our observation, we found that COL1A1 and COL1A2 expressed mainly in fibroblasts and also showed varying expression levels (Fig. S27A). Specifically, they showed a higher expression level in a subset of fibroblasts enriched with myofibroblasts (p < 0.001, Gene Set Enrichment Analysis (GSEA): Fig. S27A and B). Myofibroblasts are differentiated fibroblasts, whose dysregulation may lead to idiopathic pulmonary fibrosis (IPF) including COVID-19 lung fibrosis, a lung disease that exhibits OP [69–71]. To provide a linkage between the OP regions and IPF, we examined the scRNA-seq data from lung tissue patients with IPF and nonfibrotic controls [72]. Like in the COVID-19 data, we observed a subpopulation of fibroblasts with high expressions of COL1A1 and COL1A2 (Fig. S27C) and they were also enriched in myofibroblasts (p-values<0.001, GSEA; Fig. S27D). We further examined the differentially expressed genes in IPF fibroblast cells with high COL1A1 and COL1A2 expressions *vs.* fibroblast cells with low COL1A1 and COL1A2 expressions or nonfibrotic cells, respectively (Fig. S27E and S27F). The two lists of significant differentially expressed genes (DEGs) were similar (p- values<0.001, GSEA; Fig. S27G), suggesting that the subpopulation with low COL1A1 and COL1A2 expressions could serve as normal control. As a result, we compared the fibroblasts with a high *vs.* low COL1A1 and COL1A2 expressions in the COVID-19 scRNA-seq data and obtained the DEGs (Fig. S27H). These DEGs are enriched by two DEG lists in IPF fibroblasts with high COL1A1 and COL1A2 expressions (p< 0.001, GSEA; Fig. S27I and S27J), suggesting that these cells likely had similar molecular characteristics as those in IPF. Indeed, functional enrichment analysis using the Ingenuity Pathway Analysis reported “Pulmonary Fibrosis Idiopathic Signaling Pathways” as the 2nd top-ranked pathway (Fig. S27K). Several other top enriched pathways related to extracellular matrix remodeling were also identified. Examination of the networks of these pathways identified IL6-STAT3 and TGF-β-SMAD2/3 pathways that could directly regulate the expressions of COL1A1 and COL1A2 (Fig. 7K). Multi- color indirect immunofluorescence antibody assay (IFA) staining indeed revealed that cells with a high COL1A1 expression level also had high levels of IL6 receptor-α (IL6Rα), phospho-STAT3 (p-STAT3), TGF-βR2 and phospho-SMAD2/3 (p-SMAD2/3), confirming the activation of these pathways and their potential roles in the upregulation of COL1A1 (Fig. 7L). Both IL6 and TGF-β are highly induced, and implicated in fibrosis and lung injury [73, 74] as well as upregulation of COL1A1 and other extracellular matrix proteins in COVID-19 lung tissues [75, 76]. These results demonstrated the important roles of IL6 and TGF-β in the induction of severe COVID-19 lung OP and fibrosis. ## Discussion Using *in-situ* sequencing and our spatial single-cell analysis pipeline, we have painted an atlas of cellular and molecular signatures of lung tissues from severe COVID-19 cases. The atlas describes the cellular heterogeneity, organization, and interactions associated with inflammation, damage, and immune responses in COVID- 19 tissues. It provides molecular and cellular insights into the mechanisms underlying SARS-CoV-2 infection and the pathological manifestations in COVID-19 lungs. We detected a total of 10,414,863 transcripts in five COVID-19 and one non- COVID-19 lung tissues. Spatial single-cell transcriptomics analysis at this scale for COVID-19 and other pathological conditions is still very limited. We have overcome considerable challenges and developed a pipeline of robust, scalable, interpretable computational tools and visualization methods for targeted and exploratory analyses emphasizing on spatial discovery. The first critical component of this pipeline is a transcript-based cell segmentation strategy that integrates cell nuclei in the matching DAPI staining identified by *CellPose* with spatially localized reads using *Baysor* and iteratively segments reads into cells of potentially different sizes [10, 77]. Using this strategy, we segmented 93% of total transcripts into 1,719,459 cells including 186,659 to 470,294 cells for different tissues. A customized filter was further applied to retain only a subset of high-quality cells. Cell typing of these filtered cells based on marker genes showed that SARS- CoV-2 infected all the 18 identified cell types, a result consistent with those of our previous study [9] and others [4]. However, only a small fraction of these cells (<5%) are infected. These results agree with those of genome-wide studies that detected small numbers of viral reads in blood, lung, and nasopharyngeal samples from severe COVID-19 patients [4, 78–81]. Despite the detected low infection rates, SARS-CoV-2 inflicts similar effects on the infected and uninfected cells across the tissues as shown by differential expression analysis of genes in infected or uninfected cells in COVID-19 tissues against cells in the non-COVID-19 lung tissue (Fig. 2A). Thus, indirect effects such as those induced by the immune response, inflammation, cytokines, and complement activation might play a significant role on lung pathology in COVID-19 patients. In fact, we have observed vast pulmonary microthrombi, thrombosis, and immune infiltrations in these COVID-19 lung tissues (Fig. S1). The major parenchymal cells including ACs, fibroblasts, and VECs have the highest numbers of infected cells. The COVID-19 lung tissues have reduced numbers of ACs and SMCs, but an increased number of ECs accompanied by vast immune infiltrations of innate immune cells including MMs, granulocytes, and NK cells, and adaptive immune cells including T and B cells. Despite these global and consistent impacts of SARS-CoV-2 infection on the lungs, there are extensive spatial heterogeneities in addition to varying cell type percentages across the tissues. To reveal spatial dynamics of cellular features presented in considerable spatial discontinuity due to sparse measurements, we adopted a strategy to compute local feature distributions in the neighboring region of the cell with a fixed radius or cell number. This approach assesses the spatial gene expression maps, local infection rates, cell densities, and co-expression maps of ligand- receptor pairs. To further characterize spatial patterns of an individual feature or joint pattern relationship between two features, we extensively applied the Moran’s I statistics, a method widely adopted in geoscience but scarcely in spatial transcriptomics analysis. Using these novel strategies, we identified regions with high infection rates, which are found in regions with high cell densities expressing higher levels of SARS- CoV-2 entry-related factors including ACE2, FURIN, TMPRSS2, and NRP1. Thus, SARS-CoV-2 may preferentially infect cells expressing these factors and the infected cells could also aggregate to these high-density regions following infection. Importantly, these regions are mapped to OP (tissues 1 and 2) or regions with fibrin, hyaline membrane, and edema (tissues 3, 4, and 5), and manifest cellular and molecular patterns resembling tissue damage and wound healing including infiltration of immune cells and pattern of fibrosis. We have found that OP harbors high densities of fibroblasts expressing high levels of COL1A1 and COL1A2, a phenomenon common in fibrosis and wound healing, which are unique to OP regions unseen in other regions of COVID-19 lung tissues or in the non-COVID-19 tissue. We have identified similar cell populations enriched by myofibroblasts in 10X scRNA-seq data from COVID-19 lung samples and IPF samples, suggesting that SARS-CoV-2 might preferentially infect myofibroblasts or reprogram other cell type(s) into myofibroblasts. Pathway enrichment analysis revealed the enrichment of multiple pathways related to IPF signaling and extracellular matrix remodeling including IL6-STAT3 and TGF-b-SMAD2/3 pathways. Indeed, we have detected high levels of IL6Ra, p-STAT3, TGF-bRII, and p-SMAD2/3 in cells expressing myofibroblast markers in OP regions. These results suggest that SARS-CoV-2 infection might induce aggregation of fibroblasts, and IL6 and TGF-b to promote wound healing in OP regions. To define spatial patterns associated with COVID-19 lung pathology, we performed the NCTC analysis, or the niche analysis, commonly used for probing interactions of cells such as in immune microenvironments. However, the COVID-19 tissues contain complex pathology manifestations resulting in highly heterogenous cellular organizations and thus disparate NCTC patterns. To tackle this challenge, we applied SNMF, owing to its scalability and interpretability, to NCTCs of all cells in six tissues to reveal potentially latent features underlying these seemingly less organized NCTCs. We indeed identified 7 latent signatures of spatial cell type compositions that define normal lung structures and SARS-CoV-2 infection-induced immune niches. While many existing studies have revealed the global immune landscape of SARS-CoV-2 infection [1, 2, 5], few report spatial signatures of immune niches. We showed that these niches are mostly tissue-independent, but COVID-19 tissues presented heterogeneous distributions of these niches: while tissues 1, 2, 4, and 5 contain multiple types of niches, tissue 3 has predominately a single type. We further drew the connections between these niches and different stages of DAD. Such a linkage prompted us to apply the trajectory analysis to NCTCs of these signatures. While the trajectory analysis is most commonly applied to gene expressions to predict a pseudo-progression of the cells, here, this innovative use of trajectory analysis defined the relative severities of damages among COVID-19 tissues and constructed two different pathological routes of progression in COVID-19 patients. Both routes start with enriched “ACs with MM infiltration” and “VECs with MMs infiltration” niches but Route A is marked by increased numbers of NK cells and granulocytes, likely reflecting complications of microbial infections, while Route B is characterized by increased HIHD and OP, marked by increased fibrosis. Both routes are also correlated with multiple cytokine-cytokine receptor pairs such as CD40LG-CD40 and CXCL10-CXCR3 that are likely to mediate the chemoattraction of the immune cells (Fig. S26C and S26D). By mapping these routes to individual cells, we further revealed considerable inter and intra-tissue heterogeneity in the inferred progression of COVID-19 pathology with tissue 3 and 4 associated with Route A, tissue 2 and 5 with Route B, and tissue 1 with both routes. In summary, we have presented an atlas of spatial patterns of different cellular features that characterize SARS-CoV-2 infection, its induced immune infiltrations, inflammation, and damages in severe COVID-19 lungs, in addition to changes in gene expression in multiple cell types, including immune cells and lung parenchymal cells. These results provide insights into the spatial mechanisms at both molecular and cellular levels, which characterize the development of ARDS in COVID-19 patients. Overall, this study demonstrates the power of spatial single-cell transcriptomics and enables spatial computational analyses in the study of COVID-19 lung pathology. The developed innovative methods can benefit spatial single-cell analyses for other healthy and diseased conditions. ## Methods ### Lung tissue samples COVID-19 lung tissues were collected from 5 adults with fatal SARS-CoV-2 infection by the Autopsy Service of the Department of Pathology, Molecular and Cell- Based Medicine at the Icahn School of Medicine at Mount Sinai, and non-COVID-19 lung tissue was obtained from Pitt Biospecimen Core [9]. ### Study approval Specimens obtained at autopsy do not meet the definition of a living individual per Federal Regulations 45 CFR 46.102, and as such, research using specimens obtained at autopsy does not meet the requirements for Institutional Review Board (IRB) review or oversight under the Icahn School of Medicine Program for the Protection of Human Subjects. The University of Pittsburgh IRB determined that the study is not research involving human subjects as defined by DHHS and FDA regulations and waived of ethical oversight (STUDY20050085). ### In-situ sequencing (ISS) HS Library Preparation kit (P/N 1110-02, CARTANA AB, part of 10x Genomics) was used to prepare the library according to the manufacturer’s instruction with minor modification [82]. Four μm FFPE tissue sections were baked for 1 hour at 60°C and deparaffinized by incubating in xylene twice for 7 minutes (min) each. The sections were then rehydrated by incubating in 100% ethanol (EtOH) for 5 min, followed by 70% EtOH for 5 min, and nuclease-free deionized distilled water (SH30538.02, HyClone) for 2 min. Sections were then permeabilized by incubating in citrate buffer pH 6.0 (C9999, Sigma Aldrich) for 3 hours at 95°C. Chimeric padlock probes (from 1 custom panel (Table S1) and 4 predesigned panels: hImmune I1A, hImmune I3D, hLung L1C, hLung L2E, CARTANA AB) directly targeting RNA and containing an anchor sequence as well as a gene-specific barcode were hybridized overnight at 37°C, then ligated overnight at 30°C. Quality control of the library preparation was performed by applying anchor probes which labeled by Cy5 to simultaneously detect all rolling circle amplification products from all genes in the panel. All incubations were performed in SecureSeal™ hybridization chambers (621502, Grace Biolabs). Slow Fade Antifade Mountant (S36936, ThermoFisher) was used for mounting. All samples passed the quality control and were sent to CARTANA AB (part of 10x Genomics) for *in-situ* barcode sequencing, imaging, and data processing. Briefly, the fluorescent signals for quality control were stripped. Adapter probe pool 1 and a sequencing pool containing 4 different fluorescent labels were hybridized to the *in-situ* libraries to detect gene-specific barcodes. The scanning was fulfilled by using epifluorescence microscopy, and raw data consisting of 20x magnification images from 5 fluorescent channels (DAPI, Alexa Fluor® 488, Cy3, Cy5 and Alexa Fluor® 750) and individual z-stacks, were flattened to 2D using maximum intensity projection with a Nikon Ti2 Microscope (software NIS elements) utilizing a Zyla 4.2 camera. In total, 6 sequencing cycles were achieved for full decoding of all designed genes. After image processing, which includes image stitching, background filtering, and a sub-pixel object registration algorithm, true signals were scored based on signal intensities from individual multicolor images. The results were summarized in a CSV file and gene plots were generated using MATLAB. Each ISS spot provided a unique fluorescent barcode identifying the targeted RNA of a gene marker. The experiment generated a map of gene expressions of selected genes recorded on the natural morphology of the tissue at a single-cell level. The metadata generated from the images were further analyzed to obtain data containing the reads and their 2D (X, Y) positions of each gene, and the DAPI images. TissuuMaps [83, 84] was used to visualize spatial read locations of multiple genes against DAPI or histology images. ### Hematoxylin-eosin staining (H&E) and whole slide scanning H&E was carried out with the same slides subjected to *in-situ* sequencing using Hematoxylin & Eosin Stain Kit (H-3502, Vector Laboratories) according to the manufacturer’s instructions. The slides were then scanned with VS200 Slide Scanner (Olympus). ### Indirect immunofluorescence antibody assay (IFA) IFA was carried out as previously described [9]. Primary antibodies included CoraLite® Plus 647-conjugated smooth muscle actin (1:200, Proteintech, CL647- 67735), CoraLite® Plus 488-conjugated Collagen Type I (1:100, Proteintech, CL488- 67288), IL6Rα (1:500, Proteintech, 23457-1-AP), p-STAT3 (Tyr705) (1 to 20, Invitrogen, 710093), p-SMAD2/3 (p-SMAD2-S465/467 and p-SMAD3-S423/425, 1:500, ABclonal, AP0548), and TGF-βRII (1:10, R&D Systems, AF-241-NA). Secondary antibodies included Goat anti-Rabbit IgG (H+L) Highly Cross-Adsorbed Secondary Antibody, Alexa Fluor™ Plus 555 (1:400, Invitrogen, A32732), Donkey anti-Mouse IgG (H+L) Highly Cross-Adsorbed Secondary Antibody, Alexa Fluor™ 568 (1:400, Invitrogen, A10037). ### Cell segmentation Reads in a tissue sample were segmented into corresponding cells using the Baysor cell segmentation algorithm [10]. Baysor applies a Markov Random Field model to identify spatial clustering of reads from the same cell. It can perform cell segmentation based on read coordinates alone but can also incorporate nuclear information from DAPI staining. We first utilized the DAPI-stained images of the tissues and performed the nuclei segmentation using the *CellPose* anatomical segmentation algorithm with the default setting of parameters [77]. Segmented nuclei masks with locations were provided to Baysor to define its scale and standard deviation parameters. Baysor was able to segment 71%- 91% of reads after the first run for all samples/tissues. However, we found that among the unsegmented, the so-called noisy background reads determined by Baysor, there is a considerable portion that could still be assigned to cells visually. To address this issue, we applied Baysor again to these noisy reads from the first Baysor run. We optimized the scale and standard deviation parameter for the second Baysor run to maximize reads assignment to cells. Including the second run significantly improved the percentages of the segmented reads to 88%-95% for all samples/tissues (Fig. S3, S4 and Table S2). We implemented only two runs because the remaining unsegmented reads were mostly noisy background reads. ### Normalization and filtering Segmented gene reads assigned to cells were then filtered using multiple criteria. First, we removed cells that did not contain reads from any marker genes. Second, we removed cells with total read counts <5 or less than 4 genes with reads. Of the remaining cells, read counts for each gene within a cell were tallied and gene expressions were normalized using the *scTransform* algorithm [11] to further remove biases due to technical variability. The normalized expressions of 220 genes were used for subsequent analysis. ### Cell typing We determined the cell type of each cell using an in-house curated list of marker genes for 18 cell types (Table S1). For each cell, we calculated the average expression of marker genes for each of the 18 cell types and then assigned the cell type with the highest average expression to this cell. We further annotated infected cells as those with at least one SARS-CoV-2 read. We used UMAP to visualize the expression pattern of cells associated with the identified cell types in a lower dimension by using the visualization pipeline in the Seurat package [85]. ### Image registration The DAPI images of tissues that were used as a prior for cell segmentation had translational and rotational differences with their corresponding H&E stained images. For single-cell resolution datasets, small spatial geometric changes in the available raw morphology images could have dramatic differences in downstream analysis for studying the pathology signatures using annotated regions in the H&E image. To mitigate these spatial geometric differences, we registered the DAPI and H&E images by maximizing the phase correlation between both images. Here, we selected the DAPI image as a fixed reference because we used DAPI images as references for reads mapping and cell segmentation. To register both images, the H&E stained image is transformed by changing the scale, rotation, and translation parameters and moved across the DAPI image to improve the phase correlation of both 2D images. Once a peak correlation value was found, the optimization algorithm returned the 2D geometric transformations required to warp the H&E-stained image and registered it to the corresponding DAPI image. We applied the same image registration algorithm to all six samples used in this study. Pathology annotations were carried out on the registered H&E images. ### Cell density analysis The abundance of cells relative to the tissue context could reveal biology related to cell proliferation, cell damage, as well as other pathologies related to infection. In our spatial analysis pipeline (Fig. 1B), we defined cell density based on the number of cells within a neighborhood of a fixed radius of 200 units (200 x 0.32 μm/unit = 64 um) centering a cell, covering an area of 0.013 mm2. By calculating the number of cells within the fixed neighborhoods of each cell, we summarized the global cell density patterns across the tissue. We further presented the cell density as a contour plot by fitting a 2D kernel density estimator (KDE) on the cell density feature (Fig. 3). ### Neighboring cell type composition analysis Local cell type distributions can reveal tissue structures and different COVID-19 pathology. We assessed the local cell type distribution by computing the neighboring cell type composition (NCTC) of each cell. NCTC defined a cellular neighborhood of 200 cells and counts the number of different cell types in the neighborhood to create a neighborhood count vector for each cell. This neighborhood count vector was normalized to obtain the percentages of cell types. NCTC summarized the distributions of each cell type in the vicinity of a cell and was further used to analyze the spatial correlation of different cell types in a region-of-interest (ROI), study spatial gene expression patterns of local neighborhoods, identify local hotspots of individual cell types, and more. Furthermore, the neighborhood composition of individual cell types is spatially mapped by painting the neighborhood count/percentage of the selected cell type with a continuous “rainbow” color scheme where deeper violet to deeper reds indicates a range of low-to-high neighborhood cell composition (Fig. S10-S15). ### Spatial gene expression map SSCT captures gene expression of single cells in the context of intact tissue structures. In our study, we investigated the spatial organization of cells in tissue and their gene expression by studying their spatial gene expression maps (SGM). For each gene, its SGM is defined as the sum of the expression of that gene in its 200 nearest neighboring cells. Thus, SGM represents a spatial gradient by summarizing the regions with high or low gene expression patterns. We then painted the SGM using the continuous rainbow colormap similar to the NCTC analysis, generating the spatial gene expression maps illustrated in Fig. 2B, 2C, and 4F. By comparing the SGM with available pathology annotations, pathology signatures with relevant expression patterns can be localized to specific spatial regions in the tissue. Additionally, novel pathology signatures or ROIs can be discovered using the proposed SGM analysis. ### Local infection rate analysis We used the local infection rate analysis to study the spatial distribution of SARS-Cov-2 infection. It was similar to the NCTC analysis and computed the percentage of infected cells in a fixed physical area of 1.286 mm2 or a radius of 2,000 units (640 μm) around each cell across the tissue. These spatial rates can be further visualized using the same rainbow color scheme as NCTC (Fig. 3A). ### Ligand-receptor coexpression map Ligand-receptor interactions among cells in a defined location can define spatial patterns of immune microenvironments due to SARS-CoV-2 infection. To infer the interactions, we examine the spatial co-expression between a ligand-receptor pair. Inspired by Moran’s I spatial cross-correlation, we computed the spatial coexpression between ligand *x* and receptor *y* for cell *i* as ![Formula][1] where *x**i* = 1 (*y**j* = 1) if ligand *x* (receptor *y*) expressed in cell *i* (*j*) and 0, otherwise and 𝒩*i* defines a set of neighboring cells of cell *i*, which are cells within a radius of 1,000 units (320 μm) or a fixed physical area of 0.322 mm2 around cell *i*. Then, we visualized the coexpression in a tissue using a fixed monochrome color scheme (Fig 3F). ### Moran’s I analysis of spatial patterns Spatial global autocorrelation of one variable or more variables can be summarized using Moran’s I score. However, calculating the global spatial autocorrelation using Moran’s I across the whole tissue assumes homogeneity of the studied variable and yields only one statistic that summarizes the complete spatial pattern across the tissue disregarding their difference over space. Since our SARS- CoV-2 infected tissue samples show spatial heterogeneity, in this study, we used Local Indicators of Spatial Association (LISA) [86] to evaluate the spatial autocorrelation and the statistical significance of a study variable in each location using Local Moran’s I. Consider *x**i* as, e.g., the neighboring infection rate of a cell at location *i*, the univariate spatial autocorrelation can be found as the degree of linear association between *x**i* and a weighted average of the neighboring cells *x**j*, based on a spatial weight *W**ij*between cells at location *i* and *j*. Thus, formally, Moran’s I for each location *i* is given by ![Formula][2] where *N* defines the number of cells, *x* defines the mean of all the cells, and *m*2 is given by ![Formula][3] *I**i* ranges from -1 to 1 with -1 indicating highly negative and +1 indicating highly positive spatial autocorrelation. Such formulation can be used to identify spatial regions with high or low autocorrelation. For example, in Fig. 3B, we identified regions with high and low infection rates in each tissue. Under the same formulation, bivariate spatial autocorrelation can be used to assess spatial cross-correlation of a feature *x**i*, e.g., the neighboring infection rate of a cell at location *i* and another feature *y**j*, e.g., the cell density of a neighboring cell *j* as ![Formula][4] where *I*! defines the degree of linear association between the neighboring infection rate of a cell at location *i* and a weighted average of the cell density of the neighboring cells. Bivariate spatial associations can uncover cross-correlations between any such variable *x* at location *i* and another variable *y* at neighboring locations by ignoring correlations between *x* and *y* for cases where *i* = *j*. In our study, the bivariate LISA method is used to calculate the Local Moran’s I to identify HIHD regions described in Fig. 3E. ### Sparse Non-negative matrix factorization (SNMF) of neighboring cell type compositions SNMF was applied to extract interpretable, tissue-specific, and tissue- independent signatures from NCTC vectors from all tissue samples [63]. Let ***A*** ∈ (**0,1**)*N×M* represent the NCTC matrix of *N* cell types and *M* cells and it is factored as ![Formula][5] where ***S*** ∈ (***0,1***)*N×K* is the matrix of *K* NCTC signatures with each column denoting a signature, ***H*** ∈ (0,1)*K×M* is the signature loading matrix with *m*th column denoting the contributions of each signature to the NCTC vector of cell *m*, and ***D*** ∈ ℝ*K×K* is a scaling matrix. R package *RcppML* was used for the SNMF analysis [63], which implemented an alternating least squares (ALS) algorithm to minimize the mean squared error (MSE) between ***A*** and ***SDH***. *L*1 regularization was introduced to promote a compact, sparse signature loading. Determining the best number of signatures *K* is important for uncovering meaningful signatures. While an underestimated number could miss critical signatures, an overestimated number would include many noisy signatures. To this end, we evaluated different maximum numbers of factors from 4 to 12 and examined the meaning regarding its cell type composition and spatial distributions of loadings of each signature and chose *K* = 7 for the SNMF analysis. ### Trajectory analysis The goal of trajectory analysis is to infer a progression or ‘pseudotime’ of infection severity and associated tissue damages based on NCTC patterns in all tissues. We applied the trajectory analysis to uncover the pseudotime of 1) infection severity among tissues and 2) tissue damage defined by immune microenvironments. For 1), we used loadings of Signature 1 “Normal-like Alveoli” and Signature 3 “Broad Immune Infiltration” obtained from SNMF and defined tissue samples as clusters. For 2), we used the loadings of Signature 1 and four immune microenvironment-related signatures, *i.e.*, Signature 4-7, and performed Louvain clustering to define the clusters. The R package *Slingshot* was used to infer the trajectory [64]. The loadings, clusters, and the signatures or the 2-D UMAP matrix were fed into *Slingshot* to obtain so-called ‘pseudotime’ scores for each cell. For both cases, the cluster enriched by cells from the non-infected tissue (PBC-PR) cells was set as the origin. The pseudotime scores of the cells in each trajectory can correlate with gene expression to uncover genes that correlate with tissue damage. Additionally, ligand- receptor coexpressions that correlate with each trajectory could inform potential ligand- receptor interactions that govern immune responses. To this end, we computed the Spearman correlation of the pseudotime score with the expression of each gene per cell type in each trajectory. We also computed the correlation of ligand-receptor coexpressions with the pseudotime score of each trajectory. In both cases, to properly assess *p-values*, we performed a permutation test to obtain the empirical null distribution of the correlation coefficients and chose *p-value*<0.5 as the significant level. We then plotted heatmaps to illustrate the cell-type-wise correlations of the genes in each trajectory (Fig. S26A-D). ### Publicly available scRNA-seq datasets Two scRNA-seq datasets were obtained and processed as follows. 1. COVID-19 scRNA-seq dataset [4]. The dataset (SCP1052, lung.h5ad.gz) was downloaded from the Single Cell Portal ([https://singlecell.broadinstitute.org/single_cell](https://singlecell.broadinstitute.org/single_cell)). The processed data file lung.h5ad was accessed using the Python package Scanpy (Scanpy 1.9.1). UMAP coordinates for the processed data were accessed from the file upload.scp.X_umap.coords.txt. No follow-up processing was performed on the processed *anndata* object after it was loaded into Scanpy. 2. IPF scRNA-seq dataset [72]. The dataset (GSE135893, GSE135893_ILD_annotated_fullsize.rds.gz) was downloaded from Gene Expression Omnibus (GEO). The Seurat object file was converted into the h5ad file format using the SeuratData package (Seurat 4.3.0, SeuratData 0.2.2). After conversion, the data was loaded in Scanpy (Scanpy 1.9.1) and all cells not originating from an IPF sample or control sample were removed. ### Statistics and reproducibility Reads-based cell segmentation was performed using *Baysor* (version 0.5.2). Nuclei segmentation from DAPI images was performed using *CellPose* (version 2.1.1). All other image analyses including imaging registration were performed using MATLAB (version 2022a). Unless otherwise specified, *p*-value < 0.05 was considered significant. The differential expression analysis was performed in R (version 3.6.3) using the Wilcoxon rank-sum test with Bonferroni’s correction of multiple testing as necessary. The Moran’s I analyses for spatial patterns were performed using GeoDa (version 1.20). Pearson correlation coefficients were calculated using function *pearsonr* in Python Scipy (version 1.9.3). UMAP dimension reduction and visualization were performed using function *DimPlot* in Seurat (version 4.2.0). Trajectory analyses were performed using R package *Slingshot* (version 2.6.0). ## Data Availability The data supporting findings in this study including all the original spatial single-cell transcriptome data, segmented and filtered reads of all cells, and DAPI and H&E staining images will be deposited to Zenodo. The spatial single-cell transcriptome data will also be deposited in GEO. ## Author contributions S.-J. Gao and Y. Huang conceived, designed, supervised and managed the project. W. Meng, S.R. da Silva and L.P. Chen performed the experiments. G.L. Sica and W. Meng performed the pathological examination. A. Das, Z.T. Liu, D.M. Hasib, H. Galloway and Y. Huang developed the computational pipeline and performed the analyses. A. Das, W. Meng, Z.T. Liu, D.M. Hasib, H. Galloway, S.R. da Silva, L.P. Chen, Y.F., K.P. Rivera, M. Flores, Y.-C. Chiu, Y. Huang and S.-J. Gao interpreted the data and participated in the discussions throughout the analysis. A. Paniz-Mondolfi, C. Bryce, Z. Grimes, E.M. Sodillo and C. Cordon-Cardo obtained the COVID-19 samples. Y. Huang and S.-J. Gao wrote the manuscript with input from all the authors. All the authors read, reviewed and approved the manuscript. ## Competing interests The authors declare no competing interests. ## Data availability The data supporting findings in this study including all the original spatial single- cell transcriptome data, segmented and filtered reads of all cells, and DAPI and H&E staining images will be deposited to Zenodo. The spatial single-cell transcriptome data will also be deposited in GEO. ## Code availability All codes will be available at Zenodo. ## Supplementary information The online version contains supplementary material available at: Correspondence and requests for materials should be addressed to Shou-Jiang Gao or Yufei Huang. ## Supplemental Information ### Supplemental Tables View this table: [Table S1.](http://medrxiv.org/content/early/2023/02/26/2023.02.24.23286388/T1) Table S1. Gene annotation, cell type markers, and total cells expressing individual genes and total reads of the individual genes in each lung tissue after segmentation. View this table: [Table S2.](http://medrxiv.org/content/early/2023/02/26/2023.02.24.23286388/T2) Table S2. Summary of cell segmentation results. View this table: [Table S3.](http://medrxiv.org/content/early/2023/02/26/2023.02.24.23286388/T3) Table S3. Summary of cell typing results and SARS-CoV-2 infection status. View this table: [Table S4.](http://medrxiv.org/content/early/2023/02/26/2023.02.24.23286388/T4) Table S4. Global spatial correlations between local SARS-CoV-2 infection rates and cell densities by global Moran’s I. View this table: [Table S5.](http://medrxiv.org/content/early/2023/02/26/2023.02.24.23286388/T5) Table S5. Cell composition signatures identified by sparse non-negative matrix factorization (SNMF). View this table: [Table S6.](http://medrxiv.org/content/early/2023/02/26/2023.02.24.23286388/T6) Table S6. Spearman correlations of pseudotime Trajectory A versus gene expressions. View this table: [Table S7.](http://medrxiv.org/content/early/2023/02/26/2023.02.24.23286388/T7) Table S7. Spearman correlations of pseudotime Trajectory B versus gene expressions. ### Supplementary Figures ![Supplementary Figure 1:](http://medrxiv.org/https://www.medrxiv.org/content/medrxiv/early/2023/02/26/2023.02.24.23286388/F8.medium.gif) [Supplementary Figure 1:](http://medrxiv.org/content/early/2023/02/26/2023.02.24.23286388/F8) Supplementary Figure 1: Hematoxylin and eosin (H&E) stained images of the five COVID-19 tissues (1-2C, 2-1A, 3-1A, 4-3B, 5-3B) and the non-COVID-19 tissue (PBC- PR) along with pathology annotations are illustrated. We observe various degrees of diffuse alveolar damage (DAD) in the COVID-19 tissues. Tissues 1 and 2 (1-2C and 2- 1A) had prominent organizing pneumonia (OP) while edema, hyaline membrane, and fibrin clot were prominent in tissues 3, 4, and 5 (3-1A, 4-3B, and 5-3B). ![Supplementary Figure 2:](http://medrxiv.org/https://www.medrxiv.org/content/medrxiv/early/2023/02/26/2023.02.24.23286388/F9.medium.gif) [Supplementary Figure 2:](http://medrxiv.org/content/early/2023/02/26/2023.02.24.23286388/F9) Supplementary Figure 2: 4′,6-diamidino-2-phenylindole (DAPI) staining of all tissues under study revealed the nuclei location. We used the DAPI-stained images to facilitate cell nuclei segmentation using the CellPose algorithm. ![Supplementary Figure 3:](http://medrxiv.org/https://www.medrxiv.org/content/medrxiv/early/2023/02/26/2023.02.24.23286388/F10.medium.gif) [Supplementary Figure 3:](http://medrxiv.org/content/early/2023/02/26/2023.02.24.23286388/F10) Supplementary Figure 3: A. Bar plot showing the total number of reads and the reads that are segmented to cells using the Baysor algorithm for all tissues under study. We segmented ∼89%-95% of the reads into cells. **B.** Bar plot showing the total number of segmented cells and the identified cells with a cell-marker read. A total of 1,719,459 cells were identified across the tissues. ![Supplementary Figure 4:](http://medrxiv.org/https://www.medrxiv.org/content/medrxiv/early/2023/02/26/2023.02.24.23286388/F11.medium.gif) [Supplementary Figure 4:](http://medrxiv.org/content/early/2023/02/26/2023.02.24.23286388/F11) Supplementary Figure 4: A representative image showing cell boundary polygons of the cells segmented using the Baysor algorithm. Baysor used a binary mask of the DAPI image as a prior to guide the cell segmentation. ![Supplementary Figure 5:](http://medrxiv.org/https://www.medrxiv.org/content/medrxiv/early/2023/02/26/2023.02.24.23286388/F12.medium.gif) [Supplementary Figure 5:](http://medrxiv.org/content/early/2023/02/26/2023.02.24.23286388/F12) Supplementary Figure 5: Bar plots showing the distribution of the number of transcripts (reads) that are segmented to each cell. Over 99% of the cells harbored at least 5-15 reads. ![Supplementary Figure 6:](http://medrxiv.org/https://www.medrxiv.org/content/medrxiv/early/2023/02/26/2023.02.24.23286388/F13.medium.gif) [Supplementary Figure 6:](http://medrxiv.org/content/early/2023/02/26/2023.02.24.23286388/F13) Supplementary Figure 6: Heatmaps illustrating the average expressions of marker genes in each cell type within each tissue. In all tissues, a total of 18 cell types were identified including 11 types of parenchymal cells and 7 types of immune cells. ![Supplementary Figure 7:](http://medrxiv.org/https://www.medrxiv.org/content/medrxiv/early/2023/02/26/2023.02.24.23286388/F14.medium.gif) [Supplementary Figure 7:](http://medrxiv.org/content/early/2023/02/26/2023.02.24.23286388/F14) Supplementary Figure 7: Uniform Manifold Approximation and Projection (UMAP) of the gene expressions of individual cells revealed separated clusters for major parenchymal cells including alveolar cells (ACs) and fibroblasts. By contrast, immune cells were grouped together and mixed with vascular endothelial cells (VECs) in several clusters. The poor separation of immune cell types was largely due to highly expressed Immunoglobulin kappa light chain (IGKC) and Cathepsin L (CTSL) in these cell types, which is further illustrated in Fig. S8. ![Supplementary Figure 8:](http://medrxiv.org/https://www.medrxiv.org/content/medrxiv/early/2023/02/26/2023.02.24.23286388/F15.medium.gif) [Supplementary Figure 8:](http://medrxiv.org/content/early/2023/02/26/2023.02.24.23286388/F15) Supplementary Figure 8: The UMAP of the gene expressions of all segmented and cell-typed cells across all tissues under study is illustrated in the bottom-right panel. Here, the color of the dots indicates the tissue samples. The rest of the UMAPs illustrate the expressions of different genes colored based on their expression in each cell. We observe that IGKC and CTSL are spread across many clusters, resulting in the poor separation of immune cell types which are mixed with VECs. ![Supplementary Figure 9:](http://medrxiv.org/https://www.medrxiv.org/content/medrxiv/early/2023/02/26/2023.02.24.23286388/F16.medium.gif) [Supplementary Figure 9:](http://medrxiv.org/content/early/2023/02/26/2023.02.24.23286388/F16) Supplementary Figure 9: Spatial visualization of the segmented cells painted based on their cell types along with the location of the SARS-CoV-2 reads marked with an ‘X’. ![Supplementary Figure 10:](http://medrxiv.org/https://www.medrxiv.org/content/medrxiv/early/2023/02/26/2023.02.24.23286388/F17.medium.gif) [Supplementary Figure 10:](http://medrxiv.org/content/early/2023/02/26/2023.02.24.23286388/F17) Supplementary Figure 10: The neighborhood cell type composition (NCTC) analysis for COVID-19 tissue 1-2C is illustrated based on the cell type. ![Supplementary Figure 11:](http://medrxiv.org/https://www.medrxiv.org/content/medrxiv/early/2023/02/26/2023.02.24.23286388/F18.medium.gif) [Supplementary Figure 11:](http://medrxiv.org/content/early/2023/02/26/2023.02.24.23286388/F18) Supplementary Figure 11: The neighborhood cell type composition (NCTC) analysis for COVID-19 tissue 2-1A is illustrated based on the cell type. ![Supplementary Figure 12:](http://medrxiv.org/https://www.medrxiv.org/content/medrxiv/early/2023/02/26/2023.02.24.23286388/F19.medium.gif) [Supplementary Figure 12:](http://medrxiv.org/content/early/2023/02/26/2023.02.24.23286388/F19) Supplementary Figure 12: The neighborhood cell type composition (NCTC) analysis for COVID-19 tissue 3-1A is illustrated based on the cell type. ![Supplementary Figure 13:](http://medrxiv.org/https://www.medrxiv.org/content/medrxiv/early/2023/02/26/2023.02.24.23286388/F20.medium.gif) [Supplementary Figure 13:](http://medrxiv.org/content/early/2023/02/26/2023.02.24.23286388/F20) Supplementary Figure 13: The neighborhood cell type composition (NCTC) analysis for COVID-19 tissue 4-3B is illustrated based on the cell type. ![Supplementary Figure 14:](http://medrxiv.org/https://www.medrxiv.org/content/medrxiv/early/2023/02/26/2023.02.24.23286388/F21.medium.gif) [Supplementary Figure 14:](http://medrxiv.org/content/early/2023/02/26/2023.02.24.23286388/F21) Supplementary Figure 14: The neighborhood cell type composition (NCTC) analysis for COVID-19 tissue 5-3B is illustrated based on the cell type. ![Supplementary Figure 15:](http://medrxiv.org/https://www.medrxiv.org/content/medrxiv/early/2023/02/26/2023.02.24.23286388/F22.medium.gif) [Supplementary Figure 15:](http://medrxiv.org/content/early/2023/02/26/2023.02.24.23286388/F22) Supplementary Figure 15: The neighborhood cell type composition (NCTC) analysis for non-COVID-19 tissue PBC-PR is illustrated based on the cell type. ![Supplementary Figure 16:](http://medrxiv.org/https://www.medrxiv.org/content/medrxiv/early/2023/02/26/2023.02.24.23286388/F23.medium.gif) [Supplementary Figure 16:](http://medrxiv.org/content/early/2023/02/26/2023.02.24.23286388/F23) Supplementary Figure 16: Spatial visualization of regions with blood vessels revealed the lining of VECs along the vessels together with other vessel-associated cells including ACs, fibroblasts and smooth muscle cells (SMCs). Each dot indicates a cell, and the color of each dot indicates the cell type. ![Supplementary Figure 17:](http://medrxiv.org/https://www.medrxiv.org/content/medrxiv/early/2023/02/26/2023.02.24.23286388/F24.medium.gif) [Supplementary Figure 17:](http://medrxiv.org/content/early/2023/02/26/2023.02.24.23286388/F24) Supplementary Figure 17: Spatial visualization of a representative region with bronchioles revealed the lining of epithelial cells (ECs), smooth muscle cells (SMCs) and VECs along the bronchiole together with other bronchiole-associated cells including ACs and fibroblasts. Each dot indicates a cell, and the color of each dot indicates the cell type. ![Supplementary Figure 18:](http://medrxiv.org/https://www.medrxiv.org/content/medrxiv/early/2023/02/26/2023.02.24.23286388/F25.medium.gif) [Supplementary Figure 18:](http://medrxiv.org/content/early/2023/02/26/2023.02.24.23286388/F25) Supplementary Figure 18: Spatial visualization of the expressions of the 19 dysregulated genes in different cell types and their tissue distributions in COVID-19 tissue 1-2C are illustrated. We observe significant variations in the expression patterns suggesting their complex involvements in different aspects of SARS-CoV-2 infection and COVID-19 lung pathology. The color of each dot indicates the cell type, and the size indicates their expression level. ![Supplementary Figure 19:](http://medrxiv.org/https://www.medrxiv.org/content/medrxiv/early/2023/02/26/2023.02.24.23286388/F26.medium.gif) [Supplementary Figure 19:](http://medrxiv.org/content/early/2023/02/26/2023.02.24.23286388/F26) Supplementary Figure 19: Spatial visualization of the expressions of the 19 dysregulated genes in different cell types and their tissue distributions in COVID-19 tissue 2-1A are illustrated. We observe significant variations in the expression patterns suggesting their complex involvements in different aspects of SARS-CoV-2 infection and COVID-19 lung pathology. The color of each dot indicates the cell type, and the size indicates their expression level. ![Supplementary Figure 20:](http://medrxiv.org/https://www.medrxiv.org/content/medrxiv/early/2023/02/26/2023.02.24.23286388/F27.medium.gif) [Supplementary Figure 20:](http://medrxiv.org/content/early/2023/02/26/2023.02.24.23286388/F27) Supplementary Figure 20: Spatial visualization of the expressions of the 19 dysregulated genes in different cell types and their tissue distributions in COVID-19 tissue 3-1A are illustrated. We observe significant variations in the expression patterns suggesting their complex involvements in different aspects of SARS-CoV-2 infection and COVID-19 lung pathology. The color of each dot indicates the cell type, and the size indicates their expression level. ![Supplementary Figure 21:](http://medrxiv.org/https://www.medrxiv.org/content/medrxiv/early/2023/02/26/2023.02.24.23286388/F28.medium.gif) [Supplementary Figure 21:](http://medrxiv.org/content/early/2023/02/26/2023.02.24.23286388/F28) Supplementary Figure 21: Spatial visualization of the expressions of the 19 dysregulated genes in different cell types and their tissue distributions in COVID-19 tissue 4-3B are illustrated. We observe significant variations in the expression patterns suggesting their complex involvements in different aspects of SARS-CoV-2 infection and COVID-19 lung pathology. The color of each dot indicates the cell type, and the size indicates their expression level. ![Supplementary Figure 22:](http://medrxiv.org/https://www.medrxiv.org/content/medrxiv/early/2023/02/26/2023.02.24.23286388/F29.medium.gif) [Supplementary Figure 22:](http://medrxiv.org/content/early/2023/02/26/2023.02.24.23286388/F29) Supplementary Figure 22: Spatial visualization of the expressions of the 19 dysregulated genes in different cell types and their tissue distributions in COVID-19 tissue 5-3B are illustrated. We observe significant variations in the expression patterns suggesting their complex involvements in different aspects of SARS-CoV-2 infection and COVID-19 lung pathology. The color of each dot indicates the cell type, and the size indicates their expression level. ![Supplementary Figure 23:](http://medrxiv.org/https://www.medrxiv.org/content/medrxiv/early/2023/02/26/2023.02.24.23286388/F30.medium.gif) [Supplementary Figure 23:](http://medrxiv.org/content/early/2023/02/26/2023.02.24.23286388/F30) Supplementary Figure 23: Spatial visualization of the expressions of the 19 dysregulated genes in different cell types and their tissue distributions in non-COVID-19 tissue PBC-PR are illustrated. We observe significant variations in the expression patterns suggesting their complex involvements in different aspects of SARS-CoV-2 infection and COVID-19 lung pathology. The color of each dot indicates the cell type, and the size indicates their expression level. ![Supplementary Figure 24:](http://medrxiv.org/https://www.medrxiv.org/content/medrxiv/early/2023/02/26/2023.02.24.23286388/F31.medium.gif) [Supplementary Figure 24:](http://medrxiv.org/content/early/2023/02/26/2023.02.24.23286388/F31) Supplementary Figure 24: Bubble plots of cell-type specific differential expression analysis of the high infection versus the low infection regions of the same COVID-19 tissues is illustrated. The analysis further confirmed the upregulation of SFTPA1 or SFTPA2 in ACs and COL1A2 in fibroblasts in infected tissues. Within each row, the bubble size indicates the -log10 of the corrected p-value and the color indicates the log2 fold change of the corresponding gene. ![Supplementary Figure 25:](http://medrxiv.org/https://www.medrxiv.org/content/medrxiv/early/2023/02/26/2023.02.24.23286388/F32.medium.gif) [Supplementary Figure 25:](http://medrxiv.org/content/early/2023/02/26/2023.02.24.23286388/F32) Supplementary Figure 25: Box plots of signature loadings for each identified signature generated by the SNMF decomposition of NCTC vectors are illustrated. ![Supplementary Figure 26:](http://medrxiv.org/https://www.medrxiv.org/content/medrxiv/early/2023/02/26/2023.02.24.23286388/F33.medium.gif) [Supplementary Figure 26:](http://medrxiv.org/content/early/2023/02/26/2023.02.24.23286388/F33) Supplementary Figure 26: **A.** Heatmap illustrating the Spearman correlation of the pseudotime scores of Trajectory A and the gene expression. **B.** Heatmap illustrating the Spearman correlation of the pseudotime scores of Trajectory B and the gene expression. **C.** Heatmap illustrating the Spearman correlation of the pseudotime scores of Trajectory A and the ligand-receptor co-expression values of each cell. **D.** Heatmap illustrating the Spearman correlation of the pseudotime scores of Trajectory B and the ligand-receptor co-expression values of each cell. In all plots, deeper red indicates a high positive correlation and deeper blue indicates a high negative correlation. We also apply a p-value threshold (p-values<0.5) to remove the genes with lesser correlations and the removed genes are painted white to indicate no correlation. ![Supplementary Figure 27:](http://medrxiv.org/https://www.medrxiv.org/content/medrxiv/early/2023/02/26/2023.02.24.23286388/F34.medium.gif) [Supplementary Figure 27:](http://medrxiv.org/content/early/2023/02/26/2023.02.24.23286388/F34) Supplementary Figure 27: Analysis of COVID-19 and IPF single-cell RNA-seq datasets. A. UMAP projections showing expression of COL1A1 and COL1A2, paired with fibroblast and myofibroblast cell subtypes from single-cell RNA-seq COVID-19 tissue atlas. Violin plots showing expression of COL1A1 and COL1A2 in cells from COVID-19 patients. **B.** GSEA enrichment of myofibroblasts in cells ranked by expressions of COL1A1 and COL1A2 from the single-cell COVID-19 tissue atlas with a resulting p-value < 0.001. The significant enrichment and the large number of highly ranked hits indicate that most myofibroblasts from the single-cell COVID-19 tissue atlas have high expressions of COL1A1 and COL1A2. **C.** UMAP projections showing expressions of COL1A1 and COL1A2 next to UMAP of fibroblast cell subtypes from single-cell RNA-seq pulmonary fibrosis dataset. Violin plots showing expressions of COL1A1 and COL1A2 in IPF and control patients. The UMAPs and violin plots only represent cells from IPF patients or control. **D.** GSEA enrichment of myofibroblast in cells from IPF and control donors ranked by expressions of COL1A1 and COL1A2 with a resulting p-value < 0.001. The significant enrichment and high correlation between the myofibroblast cell type and high expressions of COL1A1 and COL1A2 indicate that most myofibroblasts from the single-cell pulmonary fibrosis dataset have high expressions of COL1A1 and COL1A2. **E.** A heatmap showing expressions of significantly differentially expressed genes from the ‘IPF High’ run in cells with high and low expressions of COL1A1 and COL1A2. **F.** A heatmap showing expressions of significantly differentially expressed genes from the ‘IPF Diagnosis’ cells with high and low expressions of COL1A1 and COL1A2. **G.** GSEA enrichment of significantly differentially expressed genes from the ‘IPF Diagnosis’ run against significantly differentially expressed genes ordered by fold change from the ‘IPF High’ run with a resulting p-value < 0.001. The high enrichment of upregulated genes and downregulated genes from the ‘IPF Diagnosis’ run in the ‘IPF High’ run indicates that the two differential expression analyses are capturing similar expression patterns. **H.** A heatmap showing expression of differentially expressed genes in cells with high and low COL1A1 and COL1A2 expressions from the single-cell COVID-19 dataset. **I.** GSEA enrichment of significantly differentially expressed genes from the ‘IPF High’ run in significantly differentially expressed genes from the COVID-19 dataset with a resulting p-value < 0.001. The high ranking of the ‘IPF High’ genes indicates that similar genes co-express with COL1A1 and COL1A2 in both COIVID-19 and IPF when considering IPF and control samples. **J.** GSEA enrichment of differentially expressed genes from ‘IPF Diagnosis’ against differentially expressed genes from the ’COVID High’ run ordered by fold change with a resulting p-value < 0.001. The high ranking of the ‘IPF Diagnosis’ genes indicates that the genes that co-express with COL1A1 and COL1A2 in COVID-19 have similar co-expression patterns in IPF. **K.** Top 10 enriched pathways from IPA pathway analysis for the differentially expressed genes from the COVID-19 single-cell RNA-seq dataset. ## Acknowledgments We thank members of Drs. Shou-Jiang Gao and Yufei Huang laboratories for technical assistance and discussions. This study was supported by grants from the National Institutes of Health (CA096512 and CA124332 to S.-J. Gao), UPMC Hillman Cancer Center Startup Funds to S.-J. Gao and Y. Huang, and in part by award P30CA047904. This research was also supported in part by the University of Pittsburgh Center for Research Computing, RRID:SCR_022735, through the resources provided. Specifically, this work used the HTC cluster, which is supported by the National Institutes of Health (S10OD028483). * Received February 24, 2023. * Revision received February 24, 2023. * Accepted February 26, 2023. * © 2023, Posted by Cold Spring Harbor Laboratory This pre-print is available under a Creative Commons License (Attribution 4.0 International), CC BY 4.0, as described at [http://creativecommons.org/licenses/by/4.0/](http://creativecommons.org/licenses/by/4.0/) ## Reference 1. 1.Ren, X., et al., COVID-19 immune features revealed by a large-scale single-cell transcriptome atlas. Cell, 2021. 184(23): p. 5838. [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1016/j.cell.2021.10.023&link_type=DOI) 2. 2.Melms, J.C., et al., A molecular single-cell lung atlas of lethal COVID-19. Nature, 2021. 595(7865): p. 114-119. 3. 3.Wilk, A.J., et al., A single-cell atlas of the peripheral immune response in patients with severe COVID-19. Nat Med, 2020. 26(7): p. 1070–1076. [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1038/s41591-020-0944-y&link_type=DOI) [PubMed](http://medrxiv.org/lookup/external-ref?access_num=32514174&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2023%2F02%2F26%2F2023.02.24.23286388.atom) 4. 4.Delorey, T.M., et al., COVID-19 tissue atlases reveal SARS-CoV-2 pathology and cellular targets. Nature, 2021. 595(7865): p. 107-113. [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1038/s41586-021-03570-8&link_type=DOI) [PubMed](http://medrxiv.org/lookup/external-ref?access_num=http://www.n&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2023%2F02%2F26%2F2023.02.24.23286388.atom) 5. 5.Stephenson, E., et al., Single-cell multi-omics analysis of the immune response in COVID-19. Nat Med, 2021. 27(5): p. 904–916. [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1038/s41591-021-01329-2&link_type=DOI) [PubMed](http://medrxiv.org/lookup/external-ref?access_num=33879890&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2023%2F02%2F26%2F2023.02.24.23286388.atom) 6. 6.Erjefalt, J.S., et al., Diffuse alveolar damage patterns reflect the immunological and molecular heterogeneity in fatal COVID-19. EBioMedicine, 2022. 83: p. 104229. 7. 7.Asp, M., J. Bergenstrahle, and J. Lundeberg, Spatially Resolved Transcriptomes- Next Generation Tools for Tissue Exploration. Bioessays, 2020. 42(10): p. e1900221. [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1002/bies.201900221&link_type=DOI) [PubMed](http://medrxiv.org/lookup/external-ref?access_num=32363691&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2023%2F02%2F26%2F2023.02.24.23286388.atom) 8. 8.Larsson, L., J. Frisen, and J. Lundeberg, Spatially resolved transcriptomics adds a new dimension to genomics. Nat Methods, 2021. 18(1): p. 15–18. [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1038/s41592-020-01038-7&link_type=DOI) 9. 9. Ramos da Silva, S., et al., Broad Severe Acute Respiratory Syndrome Coronavirus 2 Cell Tropism and Immunopathology in Lung Tissues From Fatal Coronavirus Disease 2019. J Infect Dis, 2021. 223(11): p. 1842–1854. 10. 10.Petukhov, V., et al., Cell segmentation in imaging-based spatial transcriptomics. Nat Biotechnol, 2022. 40(3): p. 345–354. 11. 11.Hafemeister, C. and R. Satija, Normalization and variance stabilization of single- cell RNA-seq data using regularized negative binomial regression. Genome Biol, 2019. 20(1): p. 296. [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1186/s13059-019-1874-1&link_type=DOI) [PubMed](http://medrxiv.org/lookup/external-ref?access_num=31870423&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2023%2F02%2F26%2F2023.02.24.23286388.atom) 12. 12.Arish, M., et al., COVID-19 immunopathology: From acute diseases to chronic sequelae. J Med Virol, 2022. 95(1): p. e28122. 13. 13.Whitsett, J.A., S.E. Wert, and T.E. Weaver, Alveolar surfactant homeostasis and the pathogenesis of pulmonary disease. Annu Rev Med, 2010. 61: p. 105–19. [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1146/annurev.med.60.041807.123500&link_type=DOI) [PubMed](http://medrxiv.org/lookup/external-ref?access_num=19824815&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2023%2F02%2F26%2F2023.02.24.23286388.atom) [Web of Science](http://medrxiv.org/lookup/external-ref?access_num=000274857600008&link_type=ISI) 14. 14.Lin, Z., et al., Genetic Association of Pulmonary Surfactant Protein Genes, SFTPA1, SFTPA2, SFTPB, SFTPC, and SFTPD With Cystic Fibrosis. Front Immunol, 2018. 9: p. 2256. [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.3389/fimmu.2018.02256&link_type=DOI) [PubMed](http://medrxiv.org/lookup/external-ref?access_num=http://www.n&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2023%2F02%2F26%2F2023.02.24.23286388.atom) 15. 15.Floros, J., et al., Human Surfactant Protein SP-A1 and SP-A2 Variants Differentially Affect the Alveolar Microenvironment, Surfactant Structure, Regulation and Function of the Alveolar Macrophage, and Animal and Human Survival Under Various Conditions. Front Immunol, 2021. 12: p. 681639. 16. 16.Mikerov, A.N., et al., Inhibition of hemagglutination activity of influenza A viruses by SP-A1 and SP-A2 variants expressed in CHO cells. Med Microbiol Immunol, 2008. 197(1): p. 9–12. [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1007/s00430-007-0051-4&link_type=DOI) [PubMed](http://medrxiv.org/lookup/external-ref?access_num=17520282&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2023%2F02%2F26%2F2023.02.24.23286388.atom) 17. 17.Garcia-Laorden, M.I., et al., The role of mannose-binding lectin in pneumococcal infection. Eur Respir J, 2013. 41(1): p. 131–9. [Abstract/FREE Full Text](http://medrxiv.org/lookup/ijlink/YTozOntzOjQ6InBhdGgiO3M6MTQ6Ii9sb29rdXAvaWpsaW5rIjtzOjU6InF1ZXJ5IjthOjQ6e3M6ODoibGlua1R5cGUiO3M6NDoiQUJTVCI7czoxMToiam91cm5hbENvZGUiO3M6MzoiZXJqIjtzOjU6InJlc2lkIjtzOjg6IjQxLzEvMTMxIjtzOjQ6ImF0b20iO3M6NTA6Ii9tZWRyeGl2L2Vhcmx5LzIwMjMvMDIvMjYvMjAyMy4wMi4yNC4yMzI4NjM4OC5hdG9tIjt9czo4OiJmcmFnbWVudCI7czowOiIiO30=) 18. 18.Sawada, K., et al., Pulmonary collectins protect macrophages against pore- forming activity of Legionella pneumophila and suppress its intracellular growth. J Biol Chem, 2010. 285(11): p. 8434–43. [Abstract/FREE Full Text](http://medrxiv.org/lookup/ijlink/YTozOntzOjQ6InBhdGgiO3M6MTQ6Ii9sb29rdXAvaWpsaW5rIjtzOjU6InF1ZXJ5IjthOjQ6e3M6ODoibGlua1R5cGUiO3M6NDoiQUJTVCI7czoxMToiam91cm5hbENvZGUiO3M6MzoiamJjIjtzOjU6InJlc2lkIjtzOjExOiIyODUvMTEvODQzNCI7czo0OiJhdG9tIjtzOjUwOiIvbWVkcnhpdi9lYXJseS8yMDIzLzAyLzI2LzIwMjMuMDIuMjQuMjMyODYzODguYXRvbSI7fXM6ODoiZnJhZ21lbnQiO3M6MDoiIjt9) 19. 19.Fassan, M., et al., Multi-Design Differential Expression Profiling of COVID-19 Lung Autopsy Specimens Reveals Significantly Deregulated Inflammatory Pathways and SFTPC Impaired Transcription. Cells, 2022. 11(6): p. 1011. 20. 20.Perkins, T.N., M.L. Donnell, and T.D. Oury, The axis of the receptor for advanced glycation endproducts in asthma and allergic airway disease. Allergy, 2021. 76(5): p. 1350–1366. 21. 21.Yue, Q., et al., Receptor for Advanced Glycation End Products (RAGE): A Pivotal Hub in Immune Diseases. Molecules, 2022. 27(15): p. 4922. 22. 22.Koerich, S., et al., Receptors for Advanced Glycation End Products (RAGE): Promising Targets Aiming at the Treatment of Neurodegenerative Conditions. Curr Neuropharmacol, 2022. Sep 22. doi: 10.2174/1570159X20666220922153903.Online ahead of print 23. 23.Yamaguchi, K., et al., Association of the RAGE/RAGE-ligand axis with interstitial lung disease and its acute exacerbation. Respir Investig, 2022. 60(4): p. 531–542. 24. 24.Hofmann, M.A., et al., RAGE mediates a novel proinflammatory axis: a central cell surface receptor for S100/calgranulin polypeptides. Cell, 1999. 97(7): p. 889–901. [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1016/S0092-8674(00)80801-6&link_type=DOI) [PubMed](http://medrxiv.org/lookup/external-ref?access_num=10399917&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2023%2F02%2F26%2F2023.02.24.23286388.atom) [Web of Science](http://medrxiv.org/lookup/external-ref?access_num=000081162800010&link_type=ISI) 25. 25.Kapandji, N., et al., Importance of Lung Epithelial Injury in COVID-19-associated Acute Respiratory Distress Syndrome: Value of Plasma Soluble Receptor for Advanced Glycation End-Products. Am J Respir Crit Care Med, 2021. 204(3): p. 359–362. [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1164/rccm.202104-1070LE&link_type=DOI) [PubMed](http://medrxiv.org/lookup/external-ref?access_num=http://www.n&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2023%2F02%2F26%2F2023.02.24.23286388.atom) 26. 26.Chiappalupi, S., et al., Hyperactivated RAGE in Comorbidities as a Risk Factor for Severe COVID-19-The Role of RAGE-RAS Crosstalk. Biomolecules, 2021. 11(6): p. 876. 27. 27.Zhang, K. and S.H. Phan, Cytokines and pulmonary fibrosis. Biol Signals, 1996. 5(4): p. 232–9. [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1159/000109195&link_type=DOI) [PubMed](http://medrxiv.org/lookup/external-ref?access_num=8891199&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2023%2F02%2F26%2F2023.02.24.23286388.atom) [Web of Science](http://medrxiv.org/lookup/external-ref?access_num=A1996VJ39700006&link_type=ISI) 28. 28.Zhu, Y., et al., Calcium in Vascular Smooth Muscle Cell Elasticity and Adhesion: Novel Insights Into the Mechanism of Action. Front Physiol, 2019. 10: p. 852. 29. 29.Ringvold, H.C. and R.A. Khalil, Protein Kinase C as Regulator of Vascular Smooth Muscle Function and Potential Target in Vascular Disorders. Adv Pharmacol, 2017. 78: p. 203–301. 30. 30.Morita, K., et al., Endothelial claudin: claudin-5/TMVCF constitutes tight junction strands in endothelial cells. J Cell Biol, 1999. 147(1): p. 185–94. [Abstract/FREE Full Text](http://medrxiv.org/lookup/ijlink/YTozOntzOjQ6InBhdGgiO3M6MTQ6Ii9sb29rdXAvaWpsaW5rIjtzOjU6InF1ZXJ5IjthOjQ6e3M6ODoibGlua1R5cGUiO3M6NDoiQUJTVCI7czoxMToiam91cm5hbENvZGUiO3M6MzoiamNiIjtzOjU6InJlc2lkIjtzOjk6IjE0Ny8xLzE4NSI7czo0OiJhdG9tIjtzOjUwOiIvbWVkcnhpdi9lYXJseS8yMDIzLzAyLzI2LzIwMjMuMDIuMjQuMjMyODYzODguYXRvbSI7fXM6ODoiZnJhZ21lbnQiO3M6MDoiIjt9) 31. 31.Schlingmann, B., S.A. Molina, and M. Koval, Claudins: Gatekeepers of lung epithelial function. Semin Cell Dev Biol, 2015. 42: p. 47–57. [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1016/j.semcdb.2015.04.009&link_type=DOI) [PubMed](http://medrxiv.org/lookup/external-ref?access_num=25951797&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2023%2F02%2F26%2F2023.02.24.23286388.atom) 32. 32.Williams, M.R., et al., Gene expression of endothelial cells due to interleukin-1 beta stimulation and neutrophil transmigration. Endothelium, 2008. 15(1): p. 73–84. [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1080/10623320802092443&link_type=DOI) [PubMed](http://medrxiv.org/lookup/external-ref?access_num=18568947&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2023%2F02%2F26%2F2023.02.24.23286388.atom) 33. 33.Hashimoto, R., et al., SARS-CoV-2 disrupts respiratory vascular barriers by suppressing Claudin-5 expression. Sci Adv, 2022. 8(38): p. eabo6783. 34. 34.Pang, L., et al., Role of caveolin-1 in human organ function and disease: friend or foe? Carcinogenesis, 2022. 43(1): p. 2–11. 35. 35.Gokani, S. and L.K. Bhatt, Caveolin-1: A Promising Therapeutic Target for Diverse Diseases. Curr Mol Pharmacol, 2022. 15(5): p. 701–715. 36. 36.Lagrange, J., et al., Alpha-2-macroglobulin in hemostasis and thrombosis: An underestimated old double-edged sword. J Thromb Haemost, 2022. 20(4): p. 806–815. 37. 37.Vandooren, J. and Y. Itoh, Alpha-2-Macroglobulin in Inflammation, Immunity and Infections. Front Immunol, 2021. 12: p. 803244. [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.3389/fimmu.2021.803244&link_type=DOI) 38. 38.Mathison, R., J.S. Davison, and A.D. Befus, Neural regulation of neutrophil involvement in pulmonary inflammation. Comp Biochem Physiol C Comp Pharmacol Toxicol, 1993. 106(1): p. 39–48. 39. 39.Sahu, S.K., et al., Emerging roles of the complement system in host-pathogen interactions. Trends Microbiol, 2022. 30(4): p. 390–402. 40. 40.Magro, C., et al., Complement associated microvascular injury and thrombosis in the pathogenesis of severe COVID-19 infection: A report of five cases. Transl Res, 2020. 220: p. 1–13. [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1016/j.trsl.2020.04.007&link_type=DOI) [PubMed](http://medrxiv.org/lookup/external-ref?access_num=http://www.n&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2023%2F02%2F26%2F2023.02.24.23286388.atom) 41. 41.Mastaglio, S., et al., The first case of COVID-19 treated with the complement C3 inhibitor AMY-101. Clin Immunol, 2020. 215: p. 108450. [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1016/j.clim.2020.108450&link_type=DOI) [PubMed](http://medrxiv.org/lookup/external-ref?access_num=32360516&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2023%2F02%2F26%2F2023.02.24.23286388.atom) 42. 42.Afzali, B., et al., The state of complement in COVID-19. Nat Rev Immunol, 2022. 22(2): p. 77–84. [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1038/s41577-021-00665-1&link_type=DOI) [PubMed](http://medrxiv.org/lookup/external-ref?access_num=http://www.n&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2023%2F02%2F26%2F2023.02.24.23286388.atom) 43. 43.Freeley, S., et al., Asparaginyl Endopeptidase (Legumain) Supports Human Th1 Induction via Cathepsin L-Mediated Intracellular C3 Activation. Front Immunol, 2018. 9: p. 2449. 44. 44.Zizzo, G., et al., Immunotherapy of COVID-19: Inside and Beyond IL-6 Signalling. Front Immunol, 2022. 13: p. 795315. 45. 45.Zhao, M.M., et al., Cathepsin L plays a key role in SARS-CoV-2 infection in humans and humanized mice and is a promising target for new drug development. Signal Transduct Target Ther, 2021. 6(1): p. 134. 46. 46.Trapani, J.A. and M.J. Smyth, Killing by cytotoxic T cells and natural killer cells: multiple granule serine proteases as initiators of DNA fragmentation. Immunol Cell Biol, 1993. 71 **(Pt** **3****)**: p. 201–8. 47. 47.Berezin, V., et al., Targeting of ECM molecules and their metabolizing enzymes and receptors for the treatment of CNS diseases. Prog Brain Res, 2014. 214: p. 353–88. [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1016/B978-0-444-63486-3.00015-3&link_type=DOI) [PubMed](http://medrxiv.org/lookup/external-ref?access_num=25410365&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2023%2F02%2F26%2F2023.02.24.23286388.atom) 48. 48.Pimentel, E., Colony-stimulating factors. Ann Clin Lab Sci, 1990. 20(1): p. 36–55. [Abstract/FREE Full Text](http://medrxiv.org/lookup/ijlink/YTozOntzOjQ6InBhdGgiO3M6MTQ6Ii9sb29rdXAvaWpsaW5rIjtzOjU6InF1ZXJ5IjthOjQ6e3M6ODoibGlua1R5cGUiO3M6NDoiQUJTVCI7czoxMToiam91cm5hbENvZGUiO3M6NDoiYWNscyI7czo1OiJyZXNpZCI7czo3OiIyMC8xLzM2IjtzOjQ6ImF0b20iO3M6NTA6Ii9tZWRyeGl2L2Vhcmx5LzIwMjMvMDIvMjYvMjAyMy4wMi4yNC4yMzI4NjM4OC5hdG9tIjt9czo4OiJmcmFnbWVudCI7czowOiIiO30=) 49. 49.Venet, F., et al., Myeloid cells in sepsis-acquired immunodeficiency. Ann N Y Acad Sci, 2021. 1499(1): p. 3–17. 50. 50.Roussel, L. and D.C. Vinh, ICOSL in host defense at epithelial barriers: lessons from ICOSLG deficiency. Curr Opin Immunol, 2021. 72: p. 21–26. 51. 51.Panneton, V., et al., Inducible T-cell co-stimulator: Signaling mechanisms in T follicular helper cells and beyond. Immunol Rev, 2019. 291(1): p. 91–103. [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1111/imr.12771&link_type=DOI) 52. 52.Mulcahy, H., et al., LST1 and NCR3 expression in autoimmune inflammation and in response to IFN-gamma, LPS and microbial infection. Immunogenetics, 2006. 57(12): p. 893–903. [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1007/s00251-005-0057-2&link_type=DOI) [PubMed](http://medrxiv.org/lookup/external-ref?access_num=16362817&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2023%2F02%2F26%2F2023.02.24.23286388.atom) [Web of Science](http://medrxiv.org/lookup/external-ref?access_num=000235127600001&link_type=ISI) 53. 53.Heidemann, J., et al., Regulated expression of leukocyte-specific transcript (LST) 1 in human intestinal inflammation. Inflamm Res, 2014. 63(7): p. 513–7. 54. 54.Chen, Z., X. Qiu, and J. Gu, Immunoglobulin expression in non-lymphoid lineage and neoplastic cells. Am J Pathol, 2009. 174(4): p. 1139–48. [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.2353/ajpath.2009.080879&link_type=DOI) [PubMed](http://medrxiv.org/lookup/external-ref?access_num=19246641&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2023%2F02%2F26%2F2023.02.24.23286388.atom) 55. 55.Vettermann, C. and M.S. Schlissel, Allelic exclusion of immunoglobulin genes: models and mechanisms. Immunol Rev, 2010. 237(1): p. 22–42. [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1111/j.1600-065X.2010.00935.x&link_type=DOI) [PubMed](http://medrxiv.org/lookup/external-ref?access_num=20727027&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2023%2F02%2F26%2F2023.02.24.23286388.atom) [Web of Science](http://medrxiv.org/lookup/external-ref?access_num=000281065300003&link_type=ISI) 56. 56.Napodano, C., et al., Mono/polyclonal free light chains as challenging biomarkers for immunological abnormalities. Adv Clin Chem, 2022. 108: p. 155–209. 57. 57.Godyna, S., M. Diaz-Ricart, and W.S. Argraves, Fibulin-1 mediates platelet adhesion via a bridge of fibrinogen. Blood, 1996. 88(7): p. 2569–77. [Abstract/FREE Full Text](http://medrxiv.org/lookup/ijlink/YTozOntzOjQ6InBhdGgiO3M6MTQ6Ii9sb29rdXAvaWpsaW5rIjtzOjU6InF1ZXJ5IjthOjQ6e3M6ODoibGlua1R5cGUiO3M6NDoiQUJTVCI7czoxMToiam91cm5hbENvZGUiO3M6MTI6ImJsb29kam91cm5hbCI7czo1OiJyZXNpZCI7czo5OiI4OC83LzI1NjkiO3M6NDoiYXRvbSI7czo1MDoiL21lZHJ4aXYvZWFybHkvMjAyMy8wMi8yNi8yMDIzLjAyLjI0LjIzMjg2Mzg4LmF0b20iO31zOjg6ImZyYWdtZW50IjtzOjA6IiI7fQ==) 58. 58.Harikrishnan, K., et al., Cell Derived Matrix Fibulin-1 Associates With Epidermal Growth Factor Receptor to Inhibit Its Activation, Localization and Function in Lung Cancer Calu-1 Cells. Front Cell Dev Biol, 2020. 8: p. 522. 59. 59.Jackson, C.B., et al., Mechanisms of SARS-CoV-2 entry into cells. Nat Rev Mol Cell Biol, 2022. 23(1): p. 3–20. [PubMed](http://medrxiv.org/lookup/external-ref?access_num=http://www.n&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2023%2F02%2F26%2F2023.02.24.23286388.atom) 60. 60.Felsenstein, S., et al., COVID-19: Immunology and treatment options. Clin Immunol, 2020. 215: p. 108448. [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1016/j.clim.2020.108448&link_type=DOI) [PubMed](http://medrxiv.org/lookup/external-ref?access_num=32353634&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2023%2F02%2F26%2F2023.02.24.23286388.atom) 61. 61.Guizani, I., et al., SARS-CoV-2 and pathological matrix remodeling mediators. Inflamm Res, 2021. 70(8): p. 847–858. 62. 62.Zhang, F., et al., SARS-CoV-2 pseudovirus infectivity and expression of viral entry-related factors ACE2, TMPRSS2, Kim-1, and NRP-1 in human cells from the respiratory, urinary, digestive, reproductive, and immune systems. J Med Virol, 2021. 93(12): p. 6671–6685. [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1002/jmv.27244&link_type=DOI) 63. 63.DeBruine, Z.J., K. Melcher, and T.J. Triche Jr, Fast and robust non-negative matrix factorization for single-cell experiments. bioRxiv, 2021: p. 2021.09. 01.458620. 64. 64.Street, K., et al., Slingshot: cell lineage and pseudotime inference for single-cell transcriptomics. BMC Genomics, 2018. 19(1): p. 477. [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1186/s12864-018-4772-0&link_type=DOI) [PubMed](http://medrxiv.org/lookup/external-ref?access_num=29914354&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2023%2F02%2F26%2F2023.02.24.23286388.atom) 65. 65.Sweeney, R.M. and D.F. McAuley, Acute respiratory distress syndrome. Lancet, 2016. 388(10058): p. 2416–2430. [PubMed](http://medrxiv.org/lookup/external-ref?access_num=http://www.n&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2023%2F02%2F26%2F2023.02.24.23286388.atom) 66. 66.Wynn, T.A., Cellular and molecular mechanisms of fibrosis. J Pathol, 2008. 214(2): p. 199–210. [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1002/path.2277&link_type=DOI) [PubMed](http://medrxiv.org/lookup/external-ref?access_num=18161745&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2023%2F02%2F26%2F2023.02.24.23286388.atom) [Web of Science](http://medrxiv.org/lookup/external-ref?access_num=000252597000008&link_type=ISI) 67. 67.Turpin, E.A., et al., Respiratory syncytial virus infection reduces lung inflammation and fibrosis in mice exposed to vanadium pentoxide. Respir Res, 2010. 11: p. 20. [PubMed](http://medrxiv.org/lookup/external-ref?access_num=20175905&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2023%2F02%2F26%2F2023.02.24.23286388.atom) 68. 68.Elewa, Y.H.A., et al., Histopathological Correlations between Mediastinal Fat- Associated Lymphoid Clusters and the Development of Lung Inflammation and Fibrosis following Bleomycin Administration in Mice. Front Immunol, 2018. 9: p. 271. 69. 69.Nagata, N., et al., Features of idiopathic pulmonary fibrosis with organizing pneumonia. Respiration, 1997. 64(5): p. 331–5. [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1159/000196700&link_type=DOI) [PubMed](http://medrxiv.org/lookup/external-ref?access_num=9311048&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2023%2F02%2F26%2F2023.02.24.23286388.atom) [Web of Science](http://medrxiv.org/lookup/external-ref?access_num=A1997XV67700003&link_type=ISI) 70. 70.Phan, S.H., The myofibroblast in pulmonary fibrosis. Chest, 2002. 122(6 Suppl): p. 286S–289S. [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1378/chest.122.6_suppl.286S&link_type=DOI) [PubMed](http://medrxiv.org/lookup/external-ref?access_num=12475801&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2023%2F02%2F26%2F2023.02.24.23286388.atom) [Web of Science](http://medrxiv.org/lookup/external-ref?access_num=000180063300006&link_type=ISI) 71. 71.Ortiz-Zapater, E., et al., Lung Fibrosis and Fibrosis in the Lungs: Is It All about Myofibroblasts? Biomedicines, 2022. 10(6). 72. 72.Habermann, A.C., et al., Single-cell RNA sequencing reveals profibrotic roles of distinct epithelial and mesenchymal lineages in pulmonary fibrosis. Sci Adv, 2020. 6(28): eaba1972. 73. 73.Laloglu, E. and H. Alay, Role of transforming growth factor-beta 1 and connective tissue growth factor levels in coronavirus disease-2019-related lung Injury: a prospective, observational, cohort study. Rev Soc Bras Med Trop, 2022. 55: p. e06152021. 74. 74. Vaz de Paula, C.B., et al., COVID-19: Immunohistochemical Analysis of TGF- beta Signaling Pathways in Pulmonary Fibrosis. Int J Mol Sci, 2021. 23(1): p. 168. 75. 75.Fan, T., et al., Discovery of 9O-Substituted Palmatine Derivatives as a New Class of anti-COL1A1 Agents via Repressing TGF-beta1/Smads and JAK1/STAT3 Pathways. Molecules, 2020. 25(4): p.773. 76. 76.Malik, M., et al., Cross-talk between Janus kinase-signal transducer and activator of transcription pathway and transforming growth factor beta pathways and increased collagen1A1 production in uterine leiomyoma cells. F S Sci, 2020. 1(2): p. 206–220. 77. 77.Stringer, C., et al., Cellpose: a generalist algorithm for cellular segmentation. Nat Methods, 2021. 18(1): p. 100–106. [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1038/s41592-020-01018-x&link_type=DOI) [PubMed](http://medrxiv.org/lookup/external-ref?access_num=33318659&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2023%2F02%2F26%2F2023.02.24.23286388.atom) 78. 78.Qi, F., et al., Single cell RNA sequencing of 13 human tissues identify cell types and receptors of human coronaviruses. Biochem Biophys Res Commun, 2020. 526(1): p. 135–140. [PubMed](http://medrxiv.org/lookup/external-ref?access_num=http://www.n&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2023%2F02%2F26%2F2023.02.24.23286388.atom) 79. 79.Sungnak, W., et al., SARS-CoV-2 entry factors are highly expressed in nasal epithelial cells together with innate immune genes. Nat Med, 2020. 26(5): p. 681–687. [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1038/s41591-020-0868-6&link_type=DOI) [PubMed](http://medrxiv.org/lookup/external-ref?access_num=http://www.n&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2023%2F02%2F26%2F2023.02.24.23286388.atom) 80. 80.Xu, G., et al., The differential immune responses to COVID-19 in peripheral and lung revealed by single-cell RNA sequencing. Cell Discov, 2020. 6: p. 73. 81. 81.Zou, X., et al., Single-cell RNA-seq data analysis on the receptor ACE2 expression reveals the potential risk of different human organs vulnerable to 2019-nCoV infection. Front Med, 2020. 14(2): p. 185–192. [PubMed](http://medrxiv.org/lookup/external-ref?access_num=http://www.n&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2023%2F02%2F26%2F2023.02.24.23286388.atom) 82. 82.Sounart, H., et al., Dual spatially resolved transcriptomics for SARS-CoV-2 host- pathogen colocalization studies in humans. bioRxiv, 2022: p. 2022.03. 14.484288. 83. 83.Pielawski, N., et al., TissUUmaps 3: Interactive visualization and quality assessment of large-scale spatial omics data. BioRxiv, 2022. doi: [https://doi.org/10.1101/2022.01.28.478131](https://doi.org/10.1101/2022.01.28.478131). 84. 84.Solorzano, L., G. Partel, and C. Wahlby, TissUUmaps: interactive visualization of large-scale spatial gene expression and tissue morphology data. Bioinformatics, 2020. 36(15): p. 4363–4365. [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1093/bioinformatics/btaa541&link_type=DOI) 85. 85.Satija, R., et al., Spatial reconstruction of single-cell gene expression data. Nat Biotechnol, 2015. 33(5): p. 495–502. [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=10.1038/nbt.3192&link_type=DOI) [PubMed](http://medrxiv.org/lookup/external-ref?access_num=25867923&link_type=MED&atom=%2Fmedrxiv%2Fearly%2F2023%2F02%2F26%2F2023.02.24.23286388.atom) 86. 86.Anselin, L., Local indicators of spatial association—LISA. Geographical analysis, 1995. 27(2): p. 93–115. [CrossRef](http://medrxiv.org/lookup/external-ref?access_num=https://doi.org/10.1111/j.1538-4632.1995.tb00338.x&link_type=DOI) [Web of Science](http://medrxiv.org/lookup/external-ref?access_num=A1995QU19800001&link_type=ISI) [1]: /embed/graphic-8.gif [2]: /embed/graphic-9.gif [3]: /embed/graphic-10.gif [4]: /embed/graphic-11.gif [5]: /embed/graphic-12.gif