香京julia种子在线播放

This study used the RNA sequencing data of normal human bronchial epithelial (NHBE) cells infected with SARS-CoV-2 (USA-WA1/2020) and mock-infected NHBE cells as control. The transcriptome data were downloaded from the Gene Expression Omnibus database (www.ncbi.nlm.nih.gov/geo/) with accession number GSE147507 (27). Each infected and controlled group has independent biological triplicate experimental data, and the reads were mapped on the human genome (hg19). The raw counts of expression data were used to identify the DEHGs between virus-infected cells compared to control with DESeq2 package version 1.22.2 in R version 3.5.2 (28). The DESeq2 takes two files as input: (a) a table of raw read counts from different samples where each row represents a gene and a column indicates the sample, and (b) an associated phenodata file describing the experimental group of each sample. DESeq2 performs internal normalization to adjust the differences in library size and library composition. For this, DESeq2 first calculates the scaling factor for each sample. Then the original read count is divided with sample scaling factor to get a normalized value. For each gene, DESeq2 uses a negative binomial distribution to model the counts and fit the normalized count data. DESeq2 uses shrinkage estimation for dispersions and fold changes and uses the Wald test to find out genes with significant differential expression between two sample groups.

A schematic diagram of our study is presented in Figure 1 . The host gene is considered as differentially expressed if the log2 Fold Change |log2FC| is > 1 and adjusted p-value is < 0.05. The expression data of the significant DEHGs we selected and transformed into Z-score (row-wise of value), then a heatmap was created with pheatmap package version 1.0.12 in R.

To understand the biological processes induced by SARS-CoV-2, the list of DEHGs were analyzed for functional enrichment, including GO of Biological Process (BP), Molecular Function (MF), and Cellular Component (CC); KEGG and BIOCARTA pathways; and DISEASE using DAVID 6.8 (https://david.ncifcrf.gov/home.jsp). The significantly enriched terms were considered when at least five genes were present, and the Benjamini adjusted p-value is less than 0.05 (p-adjust <0.05). Then the enrichment was visualized using dotplot with the ggplot2 in R.

The list of DEHGs, along with two additional host genes (ACE2 and TMPRSS2 involved in SARS-CoV-2 infection), was mapped to the STRING version 11 application (29). Although ACE2 and TMPRSS2 were not in the list of DEHGs, these two host proteins are nevertheless crucial for SARS-CoV-2 infection. Therefore, ACE2 and TMPRSS2 were taken into account during network construction to determine whether these proteins could interact with other proteins of DEHGs and influence the regulatory network. The human PPI data with medium confidence interaction score >0.400 was downloaded from STRING, and then a network was constructed and visualized in Cytoscape software version 3.7. In this network, each protein is represented by a node, while the interaction between two proteins is represented by an edge (30). Furthermore, the gene expression level (log2FC score) from DEHGs was integrated into the protein node to construct a PPI network called SiHgrn. After that, the topological structure of the SiHgrn was analyzed using Cytoscape plug-in “NetworkAnalyzer.”

The highly connected local sub-networks modules in the SiHgrn were extracted using the Cytoscape plug-in MCODE clustering algorithm (31). The functional enrichment of the highest score MCODE (Cluster 1) was analyzed for GO, KEGG, and BIOCARTA pathways, and DISEASE with DAVID 6.8. Besides, the gene set of Cluster 1 was analyzed with Molecular Signatures Database (MSigDB version 7.1) (https://www.gsea-msigdb.org/gsea/msigdb/index.jsp), which is the most famous repository of biological gene sets for performing enrichment analysis (32). Each gene set contains a collection of genes that shared a specific-biological property. In this study, MSigBD was used to identify the common underlying biological process and pathways with “hallmark gene sets” and the “positional gene set” at default parameters. The “hallmark gene sets” is the collection of 50 gene sets in which each set contains a well-defined biological state or process with coordinate expression (33). The “positional gene sets” is the collection of 299 gene sets by chromosome and cytogenetic band, and thus, it can help identify effects related to the chromosomal position. Furthermore, MSigDB was also used to determine the gene family according to their protein homology or biochemical activity.

To find the upstream TFs regulators of genes of Cluster 1, Cytoscape plug-in iRegulone (version 1.3) was used at default parameters (34). The output result was further analyzed using the in house Perl program, and a matrix table was generated where each row indicates a TF, while a column indicates a target gene. The matrix was then visualized with a heatmap using the Morpheus tool (https://software.broadinstitute.org/morpheus/). The proteins of Cluster 1 and upstream TFs are considered as a “host response signature network”.

The gene set of “host response signature network (Cluster 1 and its upstream regulators)” were analyzed to identify the biological pathways using the “Relevant Pathways” module of Network Data Exchange (NDEx version 2.4.5, 18 June 2020) (https://ndexbio.org/) (35, 36). NDEx is an open-source resource for the collection and analysis of biological network knowledge. The output result of “Relevant Pathways” analysis gave a list of biological pathways enriched for query genes. The result was sorted according to the highest similarity score between our query set and the network’s genes. Then, the top nine significantly enriched pathways having more than 5 overlapping genes were collected for further analysis. If required, figures from WikiPathways were modified with PathVisio version 3.3.0.

In order to understand how the host gene expression altered in response to SARS-CoV-2 infection, the raw read counts of host transcriptome data were preprocessed and analyzed to identify the DEHGs in SARS-CoV-2–infected cell compared to control cell with DESeq2. The boxplots showed the raw read counts, and normalized read counts across samples ( Figure S1 ). A global view of differentially expressed genes in RNA-seq datasets was visualized with the MA-plot, where M indicates the log₂ fold change on the y-axis, while A indicates the mean of the normalized counts on the x-axis ( Figure S2 ). A total of 108 DEHGs including 93 up-regulated and 15 down-regulated genes were identified with |log2FC| > 1 and adjusted p-value < 0.05. Among them, the top 15 genes showing up-regulation are CSF3, SPRR2E, CXCL5, CCL20, SPRR2D, CSF2, IL6, IFI27, IL36G, XAF1, MX1, PDZK1IP1, SAA2, IFI6, and CXCL8; while 15 genes showing down-regulation are RBM20, IFITM10, NANOS1, STON1, CXCL14, MYLK, MAP7D2, NID1, ZNF488, ATG9B, VTCN1, CYP4F3, MIR221, PPARGC1A, and METTL7A ( Table 1 ). The detailed information of up- and down-regulated genes is provided in Tables S1 and S2 , respectively. Furthermore, a heatmap of gene expression level was generated, which shows a distinct pattern of up-regulated and down-regulated genes in SARS-CoV-2–infected cells compared with mock across all biological samples ( Figure S3 ).

To understand the biological processes and key pathways altered in SARS-CoV-2–infected cell, the function and pathway enrichment analysis of the DEHGs were performed. The BP of GO analysis revealed that the up-regulated genes are primarily associated with inflammatory response, immune response, apoptotic process, type 1 interferon (IFN) signaling pathway, innate immune response, defense response to viruses, and positive regulation of NFkB transcription factor activity ( Figure 2A ). The MF of GO analysis revealed that the up-regulated genes are primarily associated with cytokine activity, growth factor activity, and chemokine activity ( Figure 2A ); while the CC of GO analysis revealed that the up-regulated genes are significantly related to extracellular space, extracellular region, and cell surface ( Figure 2A ). The KEGG pathway analysis showed the up-regulated genes are significantly enriched in TNF signaling pathway, cytokine-cytokine receptor interaction, measles, influenza A, Herpes simplex infection ( Figure 2A ). The BIOCARTA pathways showed that the up-regulated genes are significantly enriched in signal transduction through IL1R, cytokines and inflammatory response, and cells and molecules involved in the local acute inflammatory response ( Figure S4A ). The DISEASE annotation showed that the up-regulated genes are significantly enriched in type 2 diabetes, ovarian cancer, atherosclerosis, and respiratory syncytial virus bronchiolitis ( Figure S4B ). The BP, MF, and CC of GO, KEGG, BIOCARTA, and DISEASE analysis of down-regulated genes did not show any functional enrichment term. The complete results of GO, KEGG, BIOCARTA, and DISEASE annotation analyses of up-regulated host genes are available in Table S3 .

We mapped the 110 genes (108 DEHGs; ACE2 and TMPRSS2 of host genes) to the STRING database and retrieved the PPI network having 88 nodes and 563 edges. Then, the PPI network was visualized in Cytoscape software, in which each node represents a protein, while an edge represents an interaction between proteins. Furthermore, the log2FC expression level of DEHGs was integrated into the network, where the red node indicates an up-regulated gene, while the green node indicates a down-regulated gene in SARS-CoV-2–infected cell compared to control and the network is called as SiHgrn ( Figure 2B ). In the SiHgrn network, 79 nodes are up-regulated and 6 nodes are down-regulated while 3 nodes are not having gene expression level (shown in blue, identified through PPI interaction and not in the list of our DEHGs).

Though ACE2 and TMPRSS2 showed no differential expression in the SARS-CoV-2–infected cells compared to control. Interestingly, the SiHgrn network revealed that ACE2 interacts with three proteins: IL6 [log2FC=2.92], EDN1 [log2FC=1.06], and TMPRSS2; while TMPRSS2 interacts with two proteins: MX1 [log2FC=2.51] and ACE2. The structural topological of the SiHgrn network, including betweenness centrality, closeness centrality, clustering coefficient, and degree, was analyzed with Cytoscape plug-in “NetworkAnalyzer”. The complete result of network topology is given in Table S4 . The degree of node connectivity is used to determine the node size in the SiHgrn network.

Detailed analysis of the SiHgrn network with Cytoscape plug-in MCODE identified four local sub-networks ( Figure 2C ; Table 2 ). Among them, only the highest MCODE score, 16.7, Cluster 1 containing 31 genes, was selected for further study. Highly connected hub genes could influence the topology of a network and biological function. It was found that IL6 and TNF nodes have the highest degree of connectivity [degree=50] in the SiHgrn network. Other top-five hub nodes with degree of connectivity are CXCL8 [44], IL1B [40], MMP9 [35], TLR2 [34], and ICAM1 [32] ( Table S4 ). Interestingly, these all top hub genes are highly interconnected, and belong to MCODE Cluster 1 ( Figure 2C ; Table 2 ).

Function enrichment analysis of Cluster 1 genes with the BP of GO showed that they are primarily associated with immune response, inflammatory response, type 1 IFN signaling pathway, and defense response to viruses ( Figure 3A ). The MF of GO analysis revealed that Cluster 1 genes are primarily associated with chemokine activity, cytokine activity, and CXCR chemokine receptor binding ( Figure 3A ). The CC of GO of Cluster 1 genes are significantly related to extracellular space and extracellular region ( Figure 3A ). The KEGG pathway analysis showed that Cluster 1 is significantly enriched in cytokine-cytokine receptor interaction, influenza A, rheumatoid arthritis, TNF signaling pathway, measles, Herpes simplex infection ( Figure 3A ). The BIOCARTA pathways showed that the Cluster 1 genes are significantly enriched in cytokines and inflammatory response, cells and molecules involved in the local acute inflammatory response, and adhesion and diapedesis of granulocytes ( Figure S5A ). The DISEASE annotation showed that the Cluster 1 genes are significantly enriched in ovarian cancer, respiratory syncytial virus bronchiolitis, asthma, and type 2 diabetes ( Figure S5B ). The complete results of GO, KEGG, BIOCARTA, and DISEASE annotation analyses of Cluster 1 genes are available in Table S5 .

Furthermore, enrichment analysis with MSigDB hallmark gene sets showed that Cluster 1 genes are significantly enriched in several crucial pathways including, inflammatory response, genes up-regulated in response to IFN-γ and IFN-α, genes regulated by NFkB in response to TNF, genes up-regulated by IL6 via STAT3, genes up-regulated during transplant rejection and others ( Figure 3B ). While positional gene sets enrichment analysis showed that Cluster 1 genes are enriched with the location of Cytogenetic Band chr4q13 (Ensembl 99). The detailed information with the statistically significant test associated overlap matrix is presented in the Table S6 . In addition, the gene set of Cluster 1 was investigated to identify the gene family according to their protein homology or biochemical activity with MSigDB. Out of 31 genes, 15 genes (C3, CCL20, CSF2, CSF3, CXCL1, CXCL2, CXCL3, CXCL5, CXCL6, CXCL8, IL1A, IL1B, IL6, SAA1, and TNF) belong to Cytokines and growth factors; 2 genes (IRF7 and IRF9) belong to transcription factors; 4 genes (BST2, ICAM1, IFITM1, TLR2) belong to Cell differentiation markers ( Table S7 ).

Transcription factors rapidly modulate gene expression, especially in activating the host immune response against viral infection. Their identification will improve the understanding of the immunoregulatory mechanism, and a better approach to control the cytokine storm and viral defense in severe COVID-19. Therefore, we analyzed the Cluster 1 to determine its potential upstream TFs regulators using the iRegulone tool. The output result showed that most of the Cluster 1 genes, up-regulated (log2FC >1) in SARS-CoV-2 response, are under the control of TFs STAT1, STAT2, STAT3, POU2F2, and NFkB1 ( Figure 3C ).

The genes of the “host response signature network” (consist of Cluster 1 and its upstream regulators =36 query genes) were analyzed with NDEx, which returns a list of relevant pathways according to the highest similarity score. Among them, a list of top nine pathways with more than 5 overlapping genes with our query genes was selected ( Table 3 ). It was observed that the “Cytokines and Inflammatory Response” pathway significantly matched with the “host response signature network” with the highest similarity score [similarity score=0.21; the number of overlap gene=8; p-value=2.70e-10] ( Table 3 and Figure 4A ). On the other hand, “Non-genomic actions of 1,25 dihydroxyvitamin D3” pathway significantly matched with the “host response signature network” with the highest number of genes [similarity score=0.19; the number of overlap gene=11; p-value=1e-12] ( Table 3 and Figure 4B ). Figure 4B showed that TNF induced NFkB1 signaling pathway is enriched with four proteins TNF, NFkB1, CXCL8, and IL6. While, IFN-α induced Jak-STAT signaling pathway is enriched with proteins STAT1, STAT2, ILI44L, OAS1, OAS2, and OAS3. It was also observed that the gene of TLR2 is over-expressed in our study, which involves inducing the expression of CYP27B1 and VDR genes. In our study, VDR is not differentially expressed; but, CYP27B1 is up-regulated [log2FC=1.21] in the SARS-CoV-2–induced cells ( Table S1 ). Interestingly, the SiHgrn network shows interactions of CYP27B1 with IL6 and TLR2 ( Figure 2B ). Despite this, CYP27B1 is not associated with any of the MCODE clusters ( Figure 2C ). The genes of “host response signature network” are also showing enrichment with other pathways involve in immune response including “the human immune response to tuberculosis” [similarity score=0.20; the number of overlap gene=6; p-value=1.57e-8] ( Table 3 and Figure 5A ); “Photodynamic therapy-induced NF-kB survival signaling” [similarity score=0.17; the number of overlap gene=9; p-value=1e-12] ( Table 3 and Figure S6A ); “IL-10 Anti-inflammatory Signaling Pathway” [similarity score=0.16; number of overlap gene=6; p-value=4.02e-12] ( Table 3 and Figure S6B ); “Type II interferon signaling (IFN-γ)” [similarity score=0.13; the number of overlap gene=6; p-value=3.09e-7] ( Table 3 and Figure S6C ); “IL27-mediated signaling events” [similarity score=0.11; the number of overlap gene=6; p-value=3.72e-10] ( Table 3 and Figure S7A ); and “IL23-mediated signaling events” [similarity score=0.11; the number of overlap gene=7; p-value=8.21e-11] ( Table 3 and Figure S7B ). Interestingly, genes of the “host response signature network” are also showing enrichment with the “lung fibrosis” pathway [similarity score=0.12; the number of overlap gene=8; p-value=9.35e-11], which could indicate the etiology of lung damage in COVID-19 patients ( Table 3 and Figure 5B ).

To validate the identified host response pathways, a list of differentially expressed host genes was taken from a recent study conducted in 24 patients who died from COVID-19 (37). The study used the lung autopsy specimens from the 25 samples from the high viral load and 21 samples from the low viral load patients. The study identified 338 up-regulated and 5,710 down-regulated host genes in the high viral load samples compared to low viral load samples (|log2FC| >0.3 and adjusted p-value <0.01), and we referred it to as DEHGs^COVID-19 (37). Different analytical techniques, including viral gene expression, were considered for the classification of COVID-19 patients into the high and low viral load. The SARS-CoV-2 gene ORF1a is over-expressed (log2FC=9.27, adjusted p-value= 2.68e-25) in the high viral load samples compared to the low viral load samples. With the NDEx tool (28 September 2020), the up-regulated genes of DEHGs^COVID-19 identifies the “response to interferon-alpha (GO:0035455)” as a top scores relevant pathway with the similarity score=0.12; the number of overlap gene=8 (ADAR, BST2, EIF2AK2, IFIT2, IFIT3, IFITM1, IFITM3, TPR); and p-value=2.27e-6. After that, the DEHGs^COVID-19 was mapped to previously identified host response pathways with NDEx. The up-regulated and down-regulated genes with their associated pathways are given in Table 4 , and the gene expression level is provided in Table S8 .

Our analysis observed that only 4 down-regulated genes (IL1B, IL12B, TNF, and IL6) were associated with “Cytokines and Inflammatory Response” ( Table 4 and Figure 6A ). However, 10 up-regulated genes (STAT1, STAT2, IFI44L, ISG15, OAS1, OAS2, OAS3, RSAD2, PIK3R1, PIK3C2A) and 12 down-regulated genes (CAMK2B, VDR, TNF, CYP27B1, CYP24A1, CXCL8, IL6, CAMK2A, PRKCB, RXRA, PRKCA, PIK3C2B, MAPK13) were associated with the “Non-genomic actions of 1,25 dihydroxyvitamin D3” pathway ( Table 4 and Figure 6B ). The 11 up-regulated genes (STAT1, STAT2, JAK2, PSMB8, TAP1, IFI35, IFIT1, IFITM1, IFIT3, OAS1, MX1) were mapped to “the human immune response to tuberculosis” pathway ( Table 4 and Figure S8A ). Furthermore, two up-regulated genes (MT2A and EDN1) and 15 down-regulated genes (IL1B, CCR3, CXCL8, IL12B, TNF, CYSLTR2, FGF1, IL6, BMP7, MMP9, MUC5B, SFTPC, MECP2, TGFA, HMOX1) were associated with “lung fibrosis” pathway ( Table 4 and Figure S8B ). The DEHGs^COVID-19 associated with others pathways are also provided in the supplementary figures, including “Photodynamic therapy-induced NF-kB survival signaling” ( Table 4 and Figure S9A ); “IL-10 Anti-inflammatory Signaling Pathway” pathway ( Table 4 and Figure S9B ); “Type II interferon signaling (IFN-γ)” pathway ( Table 4 and Figure S9C ); “IL27-mediated signaling events” ( Table 4 and Figure S10A ); and “IL23-mediated signaling events” ( Table 4 and Figure S10B ).

Our analysis showed that the “host response signature network” is significantly associated with “cytokines and inflammatory response” pathway ( Table 3 and Figure 4A ). The gene of cytokine IL6, TNF, IL1A, IL1B, CSF2, and CSF3 were over-expressed in SARS-CoV-2 infected cells and could be responsible for proliferating hematopoietic stem cells to fight against the virus ( Figure 4A ). Other cytokine genes were also over-expressed in SARS-CoV-2 infected cells, including TNF, IL1A, and IL1B, which activate the fibroblasts and endothelial cells for inflammatory response. Figure 4A showed that CXCL2 and TNF induce neutrophil, while IL1A, IL1B, IL6, CXCL1, and TNF induce the T-helper cell. Enhance concentration of IL1A, IL1B, IL6, and TNF cause increase body temperature via the hypothalamus in COVID-19 ( Figure 4A ).

The binding of CSF2 to its receptor on the myeloid cell induces its differentiation and proliferation into monocytes and macrophages (38). However, another study showed that CSF2 required for normal pulmonary physiology, including surfactant homeostasis in the mice (39, 40). Disruption of CSF2 resulted in pulmonary alveolar proteinosis due to compromise in the development of alveolar macrophage, making lungs susceptible to infection (39, 40). During infection, CSF3 works along with IL3, IL6, and CSF2 to stimulate neutrophil granulopoiesis in the bone marrow to restore neutrophil homeostasis. The expression of CSF3 is higher in poorly controlled asthma, and inhibiting the signaling by neutralizing its receptor, CSF3R, decrease the production of mucus and hyperreactivity in the airway (41).

An effective immunological response against SARS-CoV-2 infection requires an appropriate cytokine level and an adequate number of activated T cells for reducing the viral load. A prior study found that the number of T cells was drastically decreased, while the level of TNF, IL6, and IL10 were significantly increased in severe COVID-19 patients compared to healthy control (42). Besides, the authors observed that the concentration of TNF, IL6, and IL10 was negatively correlated with the number of T cells, including total T cell, CD4+ cell, and CD8+ cells. Hence, they suggested that these cytokines play a critical role in the acute inflammatory response and lower the survival of T cells. Furthermore, the survived T cells became functionally exhausted in COVID-19, as indicated by the expression of PD-1 and Tim-3 on the T cell surface (42). The cells of healthy mice, including epithelium, endothelium, and fibroblasts, release different cytokines and involve in the immune response against pathogens (43). Our study found that several cytokine genes are over-expressed in NHBE cells infected with SARS-CoV-2. Thus, suggesting that this NHBE cell is also involved in the immunoregulatory role and is strongly implicated in the immune response against SARS-CoV-2 infection ( Figure 2A and Table S3 ), and therefore support previous findings (42, 43).

The 1,25 D is a biologically active form of vitamin D. Interestingly; prior studies observed that lower Vitamin D is associated with increased risk of autoimmune disease, inflammation, bacterial and viral infection including SARS-CoV-2 (12–15). Furthermore, an adequate amount of vitamin D could substantially reduce the risk of cytokine storm, other complications, and death in COVID-19 patients; however, its underlying mechanism is not clearly known (14, 17, 44).

Our analysis showed that the “host response signature network” is significantly associated with non-genomic actions of 1,25 D pathway with p-value 1e-12 ( Table 3 and Figure 4B ). Viral infection recognized by innate immune systems triggers two signaling cascades: (a) stimulating the production of pro-inflammatory cytokine (e.g., IL1, IL6, TNF) through NFkB1-mediated pathway; (b) stimulating the production of type I and type III IFNs through interferon regulatory factor (IRF3 and IRF7) mediated pathways (45). Based upon data analysis and literature mining, our study proposed the following two pathways through which Vitamin D could reduce the cytokine storm and enhance the antiviral response.

Only 10-20% of people infected with Mycobacterium tuberculosis have a lifetime risk of showing signs called “active” tuberculosis (TB) (58, 59). Pathogens responsible for both COVID-19 and TB show high similarity regarding its transmission mode and symptoms, including infecting the lungs, having fever, cough, and shortness of breath. A study observed that countries implemented mass vaccination programs for the Bacille Calmette-Guérin vaccine (BCG) against TB showed a significant reduction in COVID-19 mortality than those who never applied it (60). Thus, the authors suggested that BCG vaccination could protect people from COVID-19; however, the experimental evidence and underlying molecular mechanism are still lacking (60). IFN-γ is a crucial cytokine produced by CD4⁺ T cells and activates macrophage to providing resistance to TB infection (61). Low plasma concentration of IFN-γ was associated with active TB infection (62). Moreover, studies identified that a single nucleotide polymorphism (+874T/A) at the first intron of IFN-γ increases the chance to develop active TB (62, 63). A study analyzed the blood transcriptome data and identified over-expression of IFN I and IFN II inducible genes, including IFITM1, IRF1, IRF9, OAS1, MX1, STAT1, and STAT2 ( Figure 5B ), also identified up-regulated inflammatory genes in the neutrophils of active TB compared to healthy control (26). In contrast, the expression of T and B-cell specific genes were down-regulated in active tuberculosis (26, 58). In the current study, six genes (IFITM1, IRF9, MX1, OAS1, STAT1, and STAT2) of the “host response signature network” are significantly overlapping with the “human immune response to tuberculosis” pathway with p-value 1.57e-8 ( Table 3 and Figure 5B ). Furthermore, 11 genes (STAT1, STAT2, JAK2, PSMB8, TAP1, IFI35, IFIT1, IFITM1, IFIT3, OAS1, and MX1) were up-regulated in the “human immune response to tuberculosis” pathways in the high viral load compared to the low viral load COVID-19 samples ( Table 4 and Figure S8A ). Overlapping of activated genes and pathways between COVID-19 and active TB indicates both diseases are mechanistically related, which might be the reason that BCG vaccination could protect people from severe COVID-19 (60).

The complication of COVID-19 includes the development of lung fibrosis, excessive deposition of collagen and extracellular matrix (EM), destroying normal lung architecture, resulting in difficulty breathing and lung failure (64). The ARDS is considered as one of the significant factors for the lung fibrosis (64). A prior study demonstrated that SARS-CoV-1 infection promotes lung fibrosis by enhancing the TGF-β signaling and reducing the ACE2 expression (65, 66). ACE2 receptor, responsible for SARS-CoV-1 and SARS-CoV-2 infection, is primarily involved in the degradation and clearance of ANG-II. Reducing the clearance of ANG-II induces extracellular matrix deposition and lung fibrosis (66).

Our analysis showed that eight genes (CSF2, CSF3, CXCL2, CXCL8, IL1B, IL6, MMP9, and TNF) are significantly overlapping with the lung fibrosis pathway with p-value 9.35e-11 ( Table 3 and Figure 5B ). Furthermore, 15 genes (IL1B, CCR3, CXCL8, IL12B, TNF, CYSLTR2, FGF1, IL6, BMP7, MMP9, MUC5B, SFTPC, MECP2, TGFA, HMOX1) were up-regulated in the “lung fibrosis” pathways in the low viral load compared to the high viral load COVID-19 samples ( Table 4 and Figure S8A ). Thus, our study indicates the relationship between SAR-CoV-2 infection and lung fibrosis. This study also supports the previous finding that high levels of NFkB, TNF, SFTPC, MUC5B, and other proteins could promote lung fibrosis (67). Diffuse alveolar damage and lung fibrosis are frequently observed in SARS-CoV-2 patients; however, the underlying mechanism is not clearly understood (37).