and/or IL-1B) 12 leading to a "cytokine" storm coupled with the depletion of markers for platelets, natural killer cells, and dysregulation of CD 4+ and B-cell lymphocyte populations 13 . In contrast to other respiratory viral infections (e.g., ref [14] [15] [16] [17] ), SARS-CoV-2 can efficiently replicate in cells of various tissues that express angiotensin-converting enzyme 2 (ACE2) and host serine protease (TMPRSS2) 18, 19 , thus contributing to the increased transmissibility and lower-lung pathogenicity in humans 20 . This observation led us to explore mechanisms which are unique to SARS-CoV-2 infection of lung epithelial cells and cause thrombotic events.
In this study, we utilized publicly available RNA-sequencing data, GSE100457 10 , in Normal Human Bronchial Epithelial (NHBE) cell lines infected with SARS-CoV-2 (henceforth referred to as the CoV-2 dataset). In particular, we compared the CoV-2 dataset with lung epithelial cells infected with other respiratory infections, namely, respiratory syncytial virus (RSV), influenza (H1N1), and rhinovirus (RV16) (see Methods)), using a stringent study inclusion criterion.
Differentially expressed genes (DEGs) were called for each infection (with respect to their respective controls) using limma/voom, at an adjusted p-value<0.05 (see Methods). We identified a total of 339, 27, 1781 and 208 DEGs in CoV-2, RSV, H1N1 and RV16 datasets, respectively, with a more significant proportion of genes upregulated, in all cases, in response to infection (Supplementary Table 1 ). Our analysis indicated a significant over-representation of functional categories associated with immune response, across all infections (among upregulated genes) (Figure 1 ). The categories included, "response to virus", "TNFα signaling via NFkB", "interferon response" (types I/α, II/γ) and "apoptosis". However, "coagulation" and "keratinocyte differentiation pathways" were among the top pathways to be uniquely enriched within the CoV-2 dataset. Genes including MMP9, ANXA1, C1S, C3, CFB, F3, ITGB3, MAFF, MMP1, MMP9, PDGFB, PLAT, PLAU, SERPINB2, TFPI2 associated with "coagulation" and S100A7, ANXA1, TGM1, IVL, and small proline-rich proteins SPRR2A, SPRR2E, SPRR2D, associated with the "keratinocyte differentiation pathway" were among the 161 genes that were "uniquely" significantly upregulated within the CoV-2 dataset ( Figure S1 ).
We posed the question if analyzing the molecular cascades associated uniquely with SARS-CoV-2, identified above, could translate to findings that are relevant to thrombotic outcomes seen in severe/acute COVID-19 patients 2 . In order to relate the clinical and molecular phenotypes associated with SARS-CoV-2 infection, we generated a protein interaction map (strength >0.85) from DEGs identified within the CoV-2 dataset (see Methods) and annotated it with data from the comparative analysis for the three other infections (RSV, RV16 and H1N1). We subsequently clustered this network, to identify "functional modules" relevant to the pathology of SARS-CoV-2 infection. Our clustering identified nine modules with five smaller functional modules containing genes associated uniquely with SARS-CoV-2, involving angiogenesis, tubulins, keratinocytes, and small proline-rich proteins ( Figure 2 ). Three of the four remaining modules were mainly associated with immune response. These functional modules included (i)   genes associated with response to IFN type I/III such as IFITM1, IFITM3, MX1, MX2, OAS1, OAS2, IFIT1, IFIT3, and IRF9; (ii) Activation of the IL-17 pro-inflammatory cascades (IL6, IL1B, and CXCL1) via the TNFα and NFkB signaling cascades (TNF, CSF2, CSF3, NFKB1, NFKB2, RELB, and IRAK2) (iii) Signaling cascades downstream of interleukin response including JAK-STAT containing genes such as IL6R, STAT1, SOCS3, IFNGR1 and PDGFB.
Interestingly, though these modules were broadly representative of vital immune mechanisms seen in all respiratory infections 21,22 , they also highlighted certain mechanisms unique to SARS-CoV-2. For example, subversion of innate immune response is inherent to cytopathic pathogens like SARS-CoV (a related beta coronavirus), which causes an increased release of pyroptosis products (e.g. IL-1B), inducing acute inflammatory responses. The antagonism of interferon response by SARS-CoV viral proteins has been suggested to occur at multiple stages of the interferon and NFKB signaling cascades, through multiple mechanisms, including IKBKE and TRAF3 regulation, affecting downstream STAT1 associated signaling 23 . It is possible that, in the SARS-CoV-2 infection, the observed upregulation of IKBKE, TRAF3, NFKB1, NFKB2, RELB, IRAK2 at 24hpi could be a consequence of similar mechanisms 7 The fourth large module, contained genes specifically activated within the CoV-2 dataset, involved in fibrinolysis and plasminogen activation cascades (PLAT, PLAU, PLAUR, and SERPINB2), the complement activation cascade including C3, CSB, CXCL5, and platelet aggregation (platelet activation factor receptor-PTAFR (Figure 2) ). This module was particularly noteworthy given the enrichment results ( Figure 1A ) and clinical findings of COVID-19 2 . It is analyzed further in the following sections.
Plasminogen is the precursor of the serine protease plasmin, a vital component of the fibrinolytic system, essential for ensuring immune cell infiltration and cytokine production 24 . Plasminogen can be activated to plasmin by two serine peptidases, the tissue tissue-specific (T-PA) and urokinase (U-PA) plasminogen activators encoded by PLAT and PLAU, respectively. It has been previously reported for SARS-CoV and influenza infections that dysregulation of the urokinase pathway including U-PA(PLAU) and its inhibitor PAI-1(SERPINE1) might contribute to the severity of lung disease by altering the dynamics of fibrin breakdown and intra-alveolar fibrin levels and subsequently inflammation 25,26 . Here we observed a similar dynamic with SARS-CoV-2 infection albeit through activation of tissue plasminogen, PLAT, and the inhibitor SERPINB2. T-PA (which is triggered when bound to fibrin) and F3 (tissue factor, activated in CoV-2) levels are known to correlate with d-dimer levels 27 . D-dimer, a product of fibrin degradation by plasmin, is elevated in patients with COVID-19 and has been identified as a marker for disseminated intravascular coagulopathy and a worse patient prognosis 28 . These findings indicate increased fibrinolytic activity, specifically via T-PA (PLAT) activation, in SARS-CoV-2 infections.
The outcome of any viral infection is mediated through a complex interplay between viral and host proteins, which allow for a coordinated innate immune response. Plasminogen inhibitor (PAI)-1 (SERPINE1) has been reported to function as an anti-viral factor capable of inhibiting extracellular maturation of influenza particles, specifically through their action on TMPRSS2 29 .
A similar mechanism involving PAI-2 (SERPINB2) likely exists in SARS-CoV-2 infections.
Additionally, T-PA (PLAT) has been reported to interact with ORF8 protein of SARS-CoV-2 virus 7 ; however, its consequence has not yet been elucidated. We hypothesize that if this can also activated by plasmin in vivo and in vitro 30 . Additionally, degranulation of neutrophils by T-PA has been indicated as a source for increase of matrix metalloproteinase 9 (MMP-9) 31 .
Given the observed increase in MMPs, and their role in facilitating lung inflammation in ARDS (by enabling neutrophil migration and ECM breakdown) 32 , it will be essential to evaluate the impact of T-PA treatments 33 on pulmonary remodeling (via MMPs) in patients with severe/acute COVID-19.
We identified significant upregulation of several genes encoding components of the complement system, including C3, CFB and C1S, uniquely within the CoV-2 dataset, in contrast to other upper respiratory tract infections (at 24 hpi). The complement pathway is an integral part of the innate immune response and is involved in immunosurveillance for pathogen clearance (bacterial and many viral) 34 . It traditionally serves as a vital link between the innate and the adaptive immunity, driving pro-inflammatory cascades. Several of the chemokines/chemoattractants (CXCL5, CXCL6, CXCL3, and CCL20) identified within this cluster can be also be activated by signaling events precipitated by complement activation 35 . Overstimulation of these chemokines, particularly CXCL5, is known to cause destructive inflammatory lung condition in certain pathogenic models of lung disease 36 .
Plasmin also activates the complement cascade independently of established pathways (alternate, lectin, and classical pathways) by cleaving C3 and C5 to functional C3a and C5a respectively, both of which are known to be crucially involved in the inflammatory response 37 . Furthermore, there is increasing evidence for the role of complement in coagulation. C3 is shown to bind fibrinogen and fibrin with high affinity and prolonging fibrinolysis in a concentration dependent manner 38 . Additionally, studies in animal models of thrombosis have identified a plasmin driven 
We analyzed the data recently published by Blanco-Melo et al 51 to gain mechanistic insights into the pathogenicity of SARS-CoV-2. The published dataset contained multiple cell lines treated with SARS-CoV-2 including NHBE, Calu-3, A529 (with and without exogenous expression of ACE2) in addition to COVID-19 lung and normal tissue samples (available via GSE147507). Histogram of counts within normal and diseased lung samples indicated that one COVID-19 lung sample (Covid Lung 2, Figure S4 ) is a likely an outlier and was ignored from further analysis. Based on the sample clustering results of the raw counts ( Figure S2 .B), we limited our analysis to NHBE/ normal human bronchial epithelial cell lines (hence forth referred to as the CoV-2 dataset). All available series (GSE) in GEO were extracted from the gene expression omnibus with key words-"SARS-CoV", "MERS-CoV", "RSV" or "respiratory syncytial virus", "Influenza" and "Rhinovirus" for respiratory cell lines (NHBE or BEAS-2B). Since the CoV-2 transcriptional data was processed at 24 hpi (hours post infection), we chose to compare only those infections which had cell-lines at 24 hpi yielding the following series GSE3397, GSE71766, GSE100504, GSE81909, GSE27973, and GSE28904. To limit the impact of sequencing technologies, we utilized only Affymetrix, the technology with the most coverage among the series considered. Applying the above stringent inclusion criteria yielded 3 GEO series GSE3397 52 (RSV), GSE71766 53 (Influenza/H1N1 and Rhinovirus/RV16) and GSE27973 54 (Rhinovirus/RV16). We, however, did not find studies on related beta-coronaviruses in NHBE/BEAS-2B cell lines which matched our inclusion criterion.
For the sake of reproducibility, we called differentially expressed genes (DEGs) at adjusted p.value <0.05, using GEO2R for GSE3397, GSE71766 and GSE27973 comparing infected cells with their respective mocks at 24 hpi only. The DEGs called for GSE27973 were a subset of GSE71766 RV16 comparisons and subsequently ignored. Since these were microarray studies, we aggregated probes which were significantly differentially expressed, to gene names and calculated a mean fold change (Supplementary Table 1 ) and considered these for all comparisons outlined in the manuscript. We utilized the same pipeline implemented in GEO2R to call DEGs from CoV-2 data (using the limma-voom pipeline in R). Low counts were filtered using the "FilterByExpression" feature available through the edgeR package. "topTable" was used to extract all samples under our significance threshold of p.adj <0.05.We would additionally like to point out that we reanalyzed the CoV-2 data using the DESeq2 protocol as described in the original publication consistently identified similar number of DEGs as detected through limma-voom. For the comparision of COVID-19 lung biopsy sample against two healthy tissue biopsies, we utilized the exact T-test via "edgeR" to establish significance and extracted the significant genes via "topTags".
The human protein-protein interaction network (PPIN) was downloaded from STRINGdb (v 11.0) 55 for a combined edge strength of >0.85. We extracted a SARS-CoV-2 infection relevant subnetwork from this PPIN using DEGs identified in the CoV-2 dataset. The resulting CoV-2 network contained 272 nodes and 608 edges ( Figure S3 ). We additionally annotated this CoV-2 network with differential gene expression information (foldchange) identified in all infections, if present, to allow us to identify DEGs unique to CoV-2. Functionally relevant modules were extracted using the clusterMaker plugin in Cytoscape 56 . Clusters with >3 nodes were retained for further analysis, resulting in a network size of 162 nodes and 9 clusters (Figure 2 ). The clusters were named based on the most prominent terms associated with each cluster as identified using Enrichr 57 . Functional enrichment was performed using gene ontology (biological process) and mSigDB's Hallmark genesets (v7.1). All visualizations were generated via the ClusterProfiler library 58 available through R/Bioconductor. .
The entire protein-protein interaction network extracted from STRINGdb for genes differentially regulated (DEGs) in CoV-2 dataset (see Methods) is shown here. All nodes with n>3 were extracted and functionally annotated and representd in Figure 2 . This figure additionally highlights genes that were identified as DEGs within the COVID-19 lung samples (indicated with an orange node border). Red node labels indicate upregulated genes and blue node labels indicate downregulated genes. Diamond node shapes indicate DEGs identified only within CoV-2 while circle indicate DEGs identified in more than one upper respiratory tract infection. Figure 4 Analysis of previously published SARS-CoV-2 sequencing data and rationale for using NHBE cell lines and only one COVID-19 lung tissue A. Histograms of raw counts (log 2 scale, x axis) as published in GSE145708 for biopsies obtained from normal and COVID-19 lung (n=2). These results indicate likely degradation of covid lung 2 samples. We subsequently consider this sample unusable for downstream analysis of fold changes between healthy and COVID lung B. Clustering of counts across all cell lines infected with SARS-CoV-2, Influenza (and their respective controls) and one lung biopsy as obtained from the original study. Based on the clustering tree, we reasoned that the NHBE cell line better approximates lung biopsies and used them as the primary cell line of interest for this analysis