Hidden in plain sight: The effects of BCG vaccination in the COVID-19 pandemic.

To investigate the relationship between Bacille Calmette-Guérin (BCG) vaccination and SARS-CoV-2 by a bioinformatics approach, two datasets for the SARS-CoV-2 infection group and BCG-vaccinated group were downloaded. Differentially Expressed Genes were identified. Gene ontology and pathways were functionally enriched, and networking was constructed in NetworkAnalyst. Lastly, the correlation between post-BCG vaccination and COVID-19 transcriptome signatures was established. A total of 161 DEGs (113 upregulated DEGs and 48 downregulated genes) were identified in the SARS-CoV-2 group. In the pathway enrichment analysis, a cross-reference of upregulated Kyoto Encyclopedia of Genes and Genomes pathways in SARS-CoV-2 with downregulated counterparts in the BCG-vaccinated group, resulted in the intersection of 45 common pathways, accounting for 86.5% of SARS-CoV-2 upregulated pathways. Of these intersecting pathways, a vast majority were immune and inflammatory pathways with top significance in interleukin-17, tumor necrosis factor, NOD-like receptors, and nuclear factor-κB signaling pathways. Given the inverse relationship of the specific differentially expressed gene pathways highlighted in our results, the BCG-vaccine may play a protective role against COVID-19 by mounting a nonspecific immunological response and further investigation of this relationship is warranted.

INTRODUCTION

The coronavirus pandemic now affects over 2 million people with nearly 130,000 deaths reported thus far, resulting in unprecedented ramifications in global health, social infrastructure, and economic trade. Since its emergence in Wuhan, China in 2019, SARS‐coronavirus‐2 (SARS‐CoV‐2), the causative agent of the disease, has spread rapidly across the globe to over 185 countries. Curiously, the degree of intensity has varied markedly between countries—even those with similar climates, geography, and/or healthcare infrastructure. The transmission pattern does not appear to follow previous SARS‐CoV virus transmission or climatic zones based on human transmission. , , , The incidence of the COVID‐19 pandemic varies widely with the highest number of cases in the United States, Italy, China, and the lowest in regions of the world historically vulnerable to infectious diseases such as sub‐Saharan Africa. The disproportionately high morbidity and mortality of COVID‐19 infection in some countries have been linked to several factors including, but not limited to, differences in efforts to mitigate the disease spread (i.e., social distancing, travel restrictions), circulating SARS‐CoV strains, and vaccination policies.

To confront COVID‐19, it is critical to characterize the virus' biological and immunological pathways. Recent studies have uncovered the mechanism of viral entry of the SARS‐CoV‐2 relies on human angiotensin‐converting enzyme 2 (hACE2) acting as the host cell receptor. hACE2 interacts with the glycoprotein spike (S) of SARS‐CoV‐2 to form a complex structure, a receptor‐binding motif similar to some biologically related viruses, including the pathogen responsible for the Severe Acute Respiratory Syndrome (SARS). , , , , , Another study, which analyzed the transcriptional changes in the immune genes of three COVID patients, reported increased inflammatory responses to the virus, revealing increased T cell activation and cytokine expression, indicating proinflammatory pathways may be prognostic markers and/or serve as potential targets in COVID‐19 disease. Increased interleukin (IL)‐6 expression and tumor necrosis factor (TNF‐α) have also been reported, furthering the notion that cytokines are playing a key role in COVID‐induced pneumonia. Given the virus' respiratory pathology, in addition to its unique geographical distribution, some have hypothesized that protection by the Bacille Calmette‐Guérin (BCG) vaccination, a vaccine originally developed to protect against tuberculosis (TB) disease, maybe playing a critical role against COVID‐19 spread. , , , , ,

BCG is a century‐old vaccine that is given as an attenuated live strain of Mycobacterium bovis used to confer immunity against some strains of TB. Various studies have indicated that BCG vaccination has a role extending far beyond TB treatment alone, eliciting nonspecific effects (NSEs) within the innate immune system alongside the adaptive immune response. BCG vaccination has been identified for its protective role specifically against respiratory viral infections, including influenza A and respiratory syncytial virus (RSV), and is the standard therapy for certain types of bladder cancer. Although the exact mechanism remains elusive, BCG‐protection is thought to be conferred by epigenetic and immunological moderation of the immune response through the release of proinflammatory cytokines (TNF‐α, IL‐6, and interferon γ [IFN‐γ]) and the role of vitamin D. , , More specifically, one proposed mechanism dictates that during the antimicrobial response (such as in M. tb infection), toll‐like receptors (TLRs) upregulate the expression of vitamin D receptor (VDR) on immune cells and pulmonary epithelial cells. VDR then binds to calcitriol, the active form of vitamin D, and together regulates the transcriptional activity of several antimicrobial peptides such as cathelicidin and β‐defensin. , , , Additionally, vitamin D itself alters T‐cell activation and IFN‐γ stimulation, increasing the expression of several proinflammatory cytokines. Overall, the NSEs following BCG vaccination are conferred by epigenetic and transcriptional modulation of the innate immune system, as evidenced by its far‐reaching role in several viral infections beyond TB alone.

Furthermore, it has been hypothesized that countries with nationalized BCG vaccination policies show decreased morbidity and mortality to COVID‐19 when compared to those where no such uniform policy exists (such as the United States or Italy); however, these preliminary data are limited given it is yet to be peer‐reviewed and fails to account for several confounding factors such as age, testing rates, , and the accuracy of the BGC World Atlas. Nonetheless, given the safety of BCG vaccination and the well‐characterized role of NSEs, it is theorized that BCG can serve as a temporary and safe solution until a targeted vaccination becomes available. Currently, four clinical trials are already underway in Australia, the Netherlands, and the United States involving BCG versus placebo‐controlled trials in healthcare workers involved in Covid‐19 patient care (ClinicalTrial.gov; NCT04347876, NCT04327206, NCT04328441, and NCT04348370). Furthermore, there is an additional observational study in Egypt for tuberculin positivity in COVID‐19 patients versus COVID‐19 negative, BCG‐vaccinated parallel cohort (NCT04350931).

Overall, the far‐reaching effects of BCG vaccination testify to its dynamic role: Eliciting NSEs, curbing inflammation in cancer models, and reducing viremia in several distinct pathogens, including RSV, influenza A, yellow fever, and herpes simplex virus. , Altogether, this study aims to elucidate the relationship between BCG vaccination and SARS‐CoV‐2 through bioinformatic analysis of the biological and immunological pathways underlying both. Defining this relationship is the first step to understanding the role BCG vaccination may play in putative immunity against COVID‐19.

MATERIALS AND METHODS

Data source

The microarray data analyzed in this study were obtained from the Gene Expression Omnibus (https://www.ncbi.nlm.nih.gov/geo/), accession number GSE147507, published March 25, 2020. This data set analyzed the gene expression profile of a normal human bronchial epithelial (NHBE) cell line, derived from a 79‐year‐old Caucasian female, after SARS‐CoV‐2 viral infection. The following samples were analyzed; SARS‐CoV‐2 infected NHBE cells (GSM4432381, GSM4432382, and GSM4432383) compared to mock‐treated NHBE cells (GSM4432378, GSM4432379, and GSM4432380).

To compare the transcriptomic alterations in SARS‐CoV‐2 viral infection with host immune response following administration of BCG vaccine, the RNA sequencing data (GSE87186) were retrieved from the SRA database (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE87186) in the Biojupies Analysis Notebook (https://amp.pharm.mssm.edu/biojupies/analyze) and processed by ARCHS4 (all RNA‐seq and ChIP‐seq sample and signature search) pipeline (https://amp.pharm.mssm.edu/archs4). In this data set, the immunogenicity of a recombinant BCG vaccine in healthy BCG‐naïve adults and those with no prior exposure to M. tb was evaluated. We analyzed 40 samples of eight adult vaccine recipients at five different timing postvaccination; representing Days 0, 14, 28, 56, and 84 (GSM2324141–GSM2324180). A separate analysis comparing each time point to Day 0 was first performed followed by a combined analysis for all datasets (Days 14, 28, 56, and 84) compared to Day 0.

Data preprocessing and differential expression analysis

After background correction of raw expression data, mapped multiple probes to the same genes were summarized to gene levels by using the median. Genes were filtered out based on two criteria: If variance percentile rank lower than the threshold (set at 15%) or low relative abundance (average expression signal) below 5%. To ensure similar expression distributions of each sample across the entire experiment, normalization by the log2 transformation method was performed. The quality of the normalized data set was checked with the box plot and density plot, and the Benjamini and Hochberg method was selected for the multiple testing correction. Using the Limma R package available on Bioconductor (http://bioconductor.org/packages/release/bioc/html/limma.html), differentially expressed genes (DEGs) were identified in the microarray data set. Genes with a false discovery rate (FDR) < 0.05 and log2‐fold change (FC) ≥ 1.0 were considered as significantly differentially expressed and subjected to further analysis. For the RNA seq data of the BCG vaccine, gene expression signature was generated by comparing gene expression levels of the postvaccination groups with the control group on Day 0. Altered patterns of gene expression were defined through each experimental time (FC > 1, p < .05). To visually identify the similarities and differences between different cell line samples, principal component analyses (PCA) were plotted to project high‐dimensional data into lower dimensions using a linear transformation.

Functional enrichment analysis

Gene Ontology (GO) database for functional annotation of genes was used to analyze the DEGs at the functional level in terms of biological processes, molecular function, and cellular component domains. Pathway functional analysis was performed on the Kyoto Encyclopedia of Genes and Genomes (KEGG) database, Wikipathway, and Reactome pathways. The significantly overrepresented GO and pathways of the upregulated and downregulated DEGs were identified by The Search Tool for the Retrieval of Interacting Genes (STRING) database version 11.0 (https://string-db.org/) and validated in Enrichr website (http://amp.pharm.mssm.edu/Enrichr/) with p‐value < .05 as the cut‐off criterion.

Protein‐protein interaction (PPI) network construction

The STRING database was used to construct protein interaction pairs of the screened DEGs based on their function and scores. The setting was adjusted with high confidence (combined score) > 0.9 and the existence of experimental evidence. The network was visualized in Network Analyst (https://www.networkanalyst.ca/), a visual analytics platform for comprehensive gene expression profiling with the following threshold; degree distribution of 20 and betweenness at 10. The hub nodes in the PPI network were then identified based on the connectivity degree in the network statistics (number of neighbors).

Gene regulatory networks

The gene–miRNA interactome was constructed from experimentally validated miRNA–gene interaction data collected from TarBase v7.0 (http://diana.imis.athena-innovation.gr/) and miRTarBase (http://mirtarbase.mbc.nctu.edu.tw/) databases and plotted in Network Analyst. Potential transcription factors (TFs) among DEGs were first screened in the Panther (http://pantherdb.org/), JASPAR TF binding site profile database (http://jaspar.genereg.net/), and the ENCODE ChIP‐seq data (https://www.encodeproject.org/chip-seq/transcription_factor/). Additionally, upstream regulatory networks for the DEGs were built and inferred networks combining transcription factor enrichment analysis and PPI network expansion with kinase enrichment analysis were constructed in Expression2Kinase web application (http://amp.pharm.mssm.edu/X2K/).

Correlations between postvaccination host immune response and COVID transcriptome signatures

Through KEGG pathway enrichment analysis, downregulated KEGG signaling pathways following BCG vaccination were intersected with upregulated pathways in SARS‐CoV‐2 infection using Venny 2.1.0 (https://bioinfogp.cnb.csic.es/tools/venny/). Deregulated genes in each pathway were explored and the direction of expression was compared between SARS‐CoV‐2 and BCG‐vaccinated experiments in KEGG mapper (https://www.genome.jp/kegg/mapper.html).

RESULTS

Transcriptomic changes in SARS‐CoV‐2 infection

Data exploration

Raw data for the reads are provided in Table S1. After normalization, the quality was checked in density and box plots (Figure S1A,B). Mock‐treated and infected samples showed a clear discrimination pattern in PCA (Figure S1C).

DEG screening

Of 23,710 genes screened in the experiment, 7107 genes with constant values were removed. A total of 161 annotated genes were identified to be differentially expressed, which included 113 upregulated DEGs and 48 downregulated genes (Figure S1D). The top variable genes between the SARS‐CoV‐2 and mock‐treated samples are shown in Figure 1A. Using the most variable genes, heatmap (Figure 1B) showed activation of genes in the SARS‐CoV‐2 group and global downregulation in mock‐treated one. The detailed information on DEGs is listed in Table S2.

Figure 1

Functional annotations of differentially expressed genes (DEGs). (A) Three‐dimensional scatter plot for the upregulated (red) and downregulated (blue) DEGs. (B) Gene set enrichment analysis showing the top 100 variable genes between the six samples. (C) The numbers of upregulated and downregulated DEGs. The relative expression level of DEGs is stratified by the type of genes

The upregulated genes with the highest expression level include: (a) Small proline‐rich protein 2F (SPRR2F), which is an important structural protein that provides a protective barrier in stratified squamous epithelium and is associated with antimicrobial response; (b) Granulocyte‐colony stimulating factor (CSF3), a proinflammatory cytokine, has been shown to impair CD8+ T cell functionality and act as a modulator of T cell and dendritic cell functions, (c) intercellular adhesion molecule 2 (ICAM2) is expressed on bronchial epithelial, mediates adhesive interactions important for antigen‐specific immune response, NK‐cell mediated clearance, and lymphocyte recirculation, (d) S100A7A, immunogenic‐related calcium‐binding protein, regulated by TLR4, and associated with psoriasis, (e) TNFAIP3‐interacting protein 3 (TNIP3), which binds to zinc finger protein and inhibits nuclear factor‐κB (NF‐κB) activation induced by TNF‐α, TLR4, and (f) PCNA‐AS1, the antisense and regulator for proliferating cell nuclear antigen gene (Table S2).

Among the most significantly downregulated genes, the following are noted: (a) Hes Related Family BHLH Transcription Factor With YRPW Motif 2 (HEY2) is a nuclear transcription factor that represses DNA and negatively regulates miR‐146a, IL‐6, IL‐1β, and TNF‐α expression, (b) protein tyrosine phosphatase receptor type Q (PTPRQ) phosphatase which catalyzes the dephosphorylation of phosphotyrosine and phosphatidylinositol and plays roles in cellular proliferation and differentiation and is required for auditory function, and (c) taste 1 receptor member 3 (TAS1R3) encoding for G‐protein‐coupled receptor involved in taste responses (Table S2).

Functional annotations of genes

To analyze the aberrant gene expression pattern in SARS‐CoV‐2 infection, further functional analysis, and annotation for DEGs were performed (Figure 1C). Thirteen long noncoding RNAs (lncRNAs) were deregulated following the SARS‐CoV‐2 treatment. Annotation analysis revealed that PPT2‐EGFL8 readthrough (PPT2‐EGFL8, FC = −3.5169, p = 5.86E−06) was previously associated with circulating phospho‐ and sphingolipid concentrations, and others were associated with gastric, breast, and prostate cancer such as INHBA‐antisense RNA 1 (INHBA‐AS1, FC = 4.0, p = 2.41E−07), RHPN1 antisense RNA 1 (RHPN1‐AS1, FC = 3.54, p = 4.85E−06), and ST7 overlapping transcript 4 (ST7‐OT, FC = 4–4.01, p = 2.35E−07).

Genes for four miRNAs were deregulated; miR‐936 (FC = 3.88, p = 5.4E−07), miR‐23a (FC = 3.32, p = 1.8E−05), and miR‐4257 (FC = 2.79, p = 3.1E−04) were upregulated, and miR‐29b‐2 was downregulated (FC = −3.4, p = 1.12E−05). Four small nucleolar RNA (snoRNAs) genes were upregulated; including two snoRNAs, C/D box: SNORD42A (FC = 4.95, p = 1.69E−10) and SNORD12B (FC = 3.51, p = 5.86E−06), and two snoRNAs, H/ACA box: SNORA9 (FC = 2.90, p = 1.8E−04) and SNORA71C (FC = 2.63, p = 6.8E−04). In addition, 10 pseudogenes (six upregulated and four downregulated DEGs) were identified to be deregulated. Despite being nonfunctional DNA segments that resemble functional genes, some might contain inherited or acquired promoter elements and exert beneficial regulatory function (www.GeneCards.org).

As depicted in Figure 1C, DEGs included 130 protein‐coding genes (92 up and 38 down). Annotation of these genes was revealed to enclose seven transcription factors; three were activated; (a) TAL bHLH transcription factor 2 (TAL2), a basic helix‐loop‐helix transcription factor (FC = 3.74, p = 4.7E−04), (b) dachshund family transcription factor 2 (DACH2), a winged helix/forkhead transcription factor (FC = 3.68, p = .0006), and (c) AF4/FMR2 family member 2 (AFF2), a DNA‐binding transcription factor (FC = 3.40, p = .0023), and four were downregulated, (a) HEY2, basic helix‐loop‐helix transcription factor (FC = −5.24, p = 3.44E−08), and three C2H2 zinc finger transcription factor, (b) Kruppel like factor 2 (KLF2, FC = −3.85, p = 2.6E−04), (c) early growth response 2 (EGR2, FC = −2.83, p = .024), and (d) zinc finger and BTB domain containing 16 (ZBTB16, FC = −2.74, p = .0327). Furthermore, three genes for extracellular matrix structural protein were under‐expressed; (a) IgGFc‐binding protein (FCGBP, FD = −3.46, p = 8.26E−06), (b) collagen α‐3(VI) chain (COL6A3, FC = −2.93, p = .00015), and (c) complement component C1q receptor (CD93, FC = −2.74, p = .0004). Additionally, phospholipid‐transporting ATPase IM (ATP8B4), an active transporter (FC = −2.88, p = .0002) and retroviral‐like aspartic protease 1 (ASPRV1), a viral or transposable element protein, (FC = −2.83,p = .0002) were also downregulated. Two serine protease inhibitors were also downregulated; insulin‐like growth factor‐binding protein 5 (IGFBP5, FC = −2.65, p = .0006) and Serpin B10 (SERPINB10, FC = −2.81, p = .00028).

In contrast, multiple proteases were upregulated in SARS‐CoV‐2 infection compared to mock‐treated samples. Of these, SENP3‐EIF4A1 readthrough gene (SENP3‐EIF4A1, FC = 4.15, p=8.93E‐08), chymotrypsin‐like protease (CTRL‐1, FC = 2.88, p = .00019), T‐cell differentiation antigen CD6 (FC = 2.88, p = .00021), HTRA4 (FC = 2.92, p = .00016), ABHD1 (FC = 2.83, p = .0002), and cytosolic carboxypeptidase 3 AGBL3 (FC = 2.85, p = .00024). Moreover, of the upregulated DEGs, the active transporter, aquaporin (AQP7, FC = 3.07, p = .0095), the gap junction protein gamma 2 (GJC2, FC = 3.32, p = 1.87E‐05), multiple membrane traffic proteins as MX dynamin‐like GTPase 1 (MX1, FC = 2.71, p = .0004) and synaptotagmin 5 (SYT5, FC = 2.88, p = .0002), and ion channels as γ‐aminobutyric acid type A receptor subunit rho2 (GABRR2, FC = 2.74, p = .0004) and bestrophin 4 (BEST4, FC = 2.91, p = .00017) were identified (Figure 1C and Table S2).

Functional enrichment analysis of DEGs

To explore the functions of the DEGs, the 113 upregulated and 48 downregulated genes were subjected to GO and KEGG pathway enrichment analysis. As shown in Figure 2A–C, the significantly enriched GO terms for upregulated genes were mainly related to the acute inflammatory response (GO:0002526, FDR = 6.46E−08), response to virus (GO:0009615, FDR = 7.27E−08), and stress responses (GO:0006950, FDR = 8.37E−14). Among regulatory pathways, the regulation of immune response (GO:0050776, FDR = 2.23E−08) and viral genome replication (GO:0045071, FDR = 1.60E−07) were the most enriched. Lastly, cytokine activity (GO:0005125, FDR = 2.66E−09), chemokine activity (GO:0008009, FDR = 0.0018), and chemotaxis (GO:0006935, FDR = 5.07E−07) were further significantly enriched in upregulated genes (Table S3).

Figure 2

Functional enrichment analysis and gene regulatory networks. (A) Gene ontology analysis for molecular function, (B) for biological processes, (C) for cellular components, (D) The Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis. Nodes represented pathways, with color based on its significance. Four enlarged highlighted nodes represented the most significant pathways. Bipartite gene networks with these nodes showed the association with upregulated (red) and downregulated (green) genes. (E) Gene–microRNA interaction. The data source for interaction pairs: TarBase and miRTarBase. Red and green circles represented upregulated and downregulated differentially expressed genes, respectively. Blue squares represented the microRNAs. The full network is shown with extracted nodes (ellipse) for four significant KEGG pathways. (F) Transcription factor enrichment analysis showing putative transcription factors that most likely to regulate the differences in gene expression. (G) Upstream regulatory network that connects the enriched transcription factors to kinases through known protein–protein interactions. (H) Kinase enrichment analysis. Candidate enriched protein kinases that most likely regulate the formation of the identified transcriptional complexes. They are ranked based on the overlap between known kinase‐substrate phosphorylation interactions and the proteins in the protein–protein interaction subnetwork created in (G)

KEGG pathway analysis (Figure 2D) showed that upregulated genes were mainly enriched in the IL‐17 signaling pathway (KEGG ID: hsa04657, FDR = 1.39E−15), TNF signaling pathway (hsa04668, FDR = 5.34E−15), NOD‐like receptor (NLR) signaling pathway (hsa04621, FDR = 4.07E−10), cytokine–cytokine receptor interaction (hsa04060, FDR = 2.08E−07), Jak‐STAT signaling pathway (hsa04630, FDR = 5.38E−06), NF‐κB signaling pathway (hsa04064, FDR = 1.33E−07), and influenza A (hsa05164, FDR = 1.10E−07). Other viral pathways included Kaposi's sarcoma‐associated herpesvirus infection (hsa05167), herpes simplex infection (hsa05168), measles (hsa05162), hepatitis C (hsa05160), and Epstein‐Barr virus infection (hsa05169) (Table S4). In contrast, no significant GO terms or pathways were significant for the downregulated genes.

As depicted in Figure 2D, upregulated pathways are highly connected with inflammatory cytokines exhibiting significant crosstalk in particular. IL‐6 (also known as B‐stimulatory factor‐2) is involved in the final differentiation of B cells into immunoglobulin‐secreting cells and the induction of acute‐phase reactants. IL‐6 has sequence similarity with granulocyte colony‐stimulating factors 2 and 3 (CSF2 and CSF3), which are also highly expressed following SARS‐CoV‐2 infection. Furthermore, several interleukins (IL17C, IL19, IL36A, and IL36G) were noted to serve as proinflammatory cytokines for the regulation of dendritic cells and T cells, as well as interleukin‐2 receptor subunit gamma (IL2RG), which encodes a common gamma chain essential for IL‐receptor function. Overexpression of two members of the CC chemokine family was observed: C–C Motif chemokine ligand 2 (CCL2) and C–C motif chemokine ligand 2 (CCL20), whose function was to induce the migration of monocytes and lymphocytes, respectively (www.GeneCards.org).

Gene regulatory networks

The gene–miRNA interactions network was constructed (Figure 2E). It consisted of 77 seeds (significant genes), 1136 nodes, and 1776 edges mapped to the corresponding molecular interaction databases. KEGG pathway enrichment analysis of the network was significant for four pathways; namely, IL‐17 signaling pathway (p = 7E−05, Hits = 5/93), TNF signaling pathway (p = .0017, Hits = 4/110), cytokine–cytokine receptor pathway (p = .012, Hits = 5/294), and influenza A (p = .044, Hits = 3/167). After extraction of nodes relevant to these immune‐related pathways, the densely connected microRNAs in the cluster module included: miR‐26b‐5p, miR‐26a‐5p, miR‐124‐3p, miR‐7‐5p, miR‐17‐5p, miR‐335‐5p, mirR‐24‐3p, miR‐203a‐3p, and miR‐122‐5p.

Upstream regulatory cell signaling layers responsible for the observed pattern in gene expression following SARS‐CoV‐2 infection is depicted in Figure 2F–H. The network included integrated transcription factors and kinases with PPI. The top enriched transcription factors were RELA proto‐oncogene (p = 1.55E−06), BCL3 (p = .0014), TP63 (p = 2.9E−05), and VDR (p = .0051). RELA, an NF‐κB subunit, has a key role in mediating inflammation, differentiation, cell growth, tumorigenesis, and apoptosis (www.GeneCards.org). It is associated with 12 overlapping targets in the DEGs list (IER3, STAT1, NFKB2, ICAM1, DRAM1, NEDD9, IL32, NFKBIA, TNFAIP3, NFKBIZ, and TNIP1). Another NF‐κB regulator, BCL3, has a paradoxical effect depending on its subcellular localization: In the nucleus, it regulates transcriptional activation of NF‐κB target gene; however, in the cytoplasm, it inhibits the nuclear translocation of the NF‐κB p50 subunit (www.GeneCards.org). From the DEGs defined following SARS‐CoV‐2 infection, six genes (IER3, STAT1, NFKB2, NFKBIA, IRF9, and BIRC3) were targeted by BCL3. Lastly, the nuclear receptor for calcitriol, VDR, was significantly associated with the upregulated DEGs (NFKB2, NFKBIA, and IRAK2).

TB pathway

Gene set enrichment analysis (GSEA) showed an over‐representation of SARS‐CoV‐2 genes in the TB pathway (ID: hsa05152, p = 6.09E−04). In this pathway, 116 DEGs (of 179 genes) were upregulated (Figure S2). IL6, TLR9, and TNF were at the top of that list. VDR and CEPB are further noted for their critical regulation of inflammatory and immune responses, including SARS‐CoV‐2 infected cells.

Transcriptomic changes following BCG vaccine administration

DEGs

After preprocessing of the expression data from the four different stages postvaccination (Figure S3), 696, 66, 49, and 80 DEGs were detected at 14, 28, 56, and 84 days, respectively. As depicted in Figure 3A, THE principal component analysis showed the clustering of samples with discrete demarcation from the initial stage at Day 0 (controls). A heatmap of the most variable genes demonstrated an unclear difference in the transcriptomic pattern at different time points (Figure 3B). As shown in Figure 3C, the volcano plot illustrated deregulated genes. A total of 696, 66, 49, and 80 DEGs, respectively, were identified (Figure 3D). The altered expression patterns across four different time points showed host immune responses to be constantly changing post‐BCG vaccination. Stratified analysis according to the direction of DEGs showed similar results with no overlapping between up and down DEGs at each stage (Figure 3E,F).

Figure 3

Differentially expressed genes (DEGs) and functional enrichment analysis following Bacille Calmette‐Guérin (BCG) vaccination. (A) Principal component analysis after normalization, showing a cluster of samples postvaccination on Days 14 through 84. (B) Clustergram showing the hierarchical clustering of the top 2500 variable genes among the five groups. (C) Volcano plot representing log2‐fold change and −log10 (adjusted p‐value). (D) Number of DEGs at each stage postvaccination (Days 14 vs. 0, 28 vs. 0, 56 vs. 0, and 84 vs. 0) (E, F) Venn diagram showing the intersection between upregulated and downregulated DEGs at a different stage. The numbers of upregulated and downregulated DEGs are demonstrated in blue and red boxes, respectively. (G) Protein–protein interaction (PPI) network for DEGs following BCG vaccination. String interactome for PPI showing upregulated (red nodes) and downregulated (green nodes) genes. (H) A cluster of inhibited Kyoto Encyclopedia of Genes and Genomes pathways. Circled in blue the gene list of the selected pathways. (I) Top significant pathways inhibited following BCG vaccination (enriched in the cluster)

Functional enrichment analysis

Activated and inhibited GO and pathway terms for each stage are listed in detail in Tables S5–S12. GO enrichment analysis was carried out in three categories, including biological processes (BP), molecular functions (MF), and cellular components (CC). The most significant terms for BP, MF, and CC were the translation (GO:0006412, FDR = 8.59E−10), RNA binding (GO:0003723, FDR = 6.98E−08), and cytosolic ribosome (GO:0022626, FDR = 1.30E−08) at Day 14; short‐chain fatty acid catabolic process (GO:0019626, FDR = 0.0008), azole transmembrane transporter activity (GO:1901474, FDR = 0.0088), and nucleolar ribonuclease P complex (GO:0005655, FDR = 0.1409) at Day 28; protein import into peroxisome matrix (GO:0016558, FDR = 0.0029), β‐1,3‐galactosyltransferase activity (GO:0048531, FDR = 0.0046), and manchette (GO:0002177, FDR = 0.0121) at Day 56; and positive regulation of T‐helper 2 cell differentiation (GO:0045630, FDR = 0.0088), primary lysosome (GO:0005766, FDR = 0.0246), and alanine transmembrane transporter activity (GO:0022858, FDR = 0.0121) at Day 84 postvaccination.

Pathway enrichment analysis showed multiple activated pathways at different stages. The most significant upregulated pathways at an early stage (Day‐14 and Day‐28 groups) were (a) ribosome (hsa03010, FDR = 3.83E−12) including multiple mitochondrial ribosomal proteins (as MRPL17, MRPS21, RPL23A, MRPL24, MRPL34, and MRPL36), and ribosomal proteins for small subunits (RPS14, RPS15, RPS16, RPS19, and RPS28), and large subunits (RPL18A, RPL36, RPLP2, RPL13, and RPL15) and (b) propionyl‐CoA catabolism (R‐HSA‐71032, FDR = 0.0001) via modulating methyl malonyl CoA epimerase (MCEE), methylmalonic aciduria type A (MMAA), and propionyl‐CoA carboxylase (PCCB) genes. However, at later stages, the top upregulated pathways were aminoacyl‐tRNA biosynthesis (hsa00970, FDR = 0.0246), scavenging of heme from plasma (R‐HSA‐2168880, FDR = 0.0056), and heme biosynthesis (WP561, FDR = 0.020).

In contrast, of the inhibited pathways after 2 weeks of vaccination were cytokine–cytokine receptor interaction (hsa04060, FDR = 2.63E−04), JAK‐STAT signaling pathway (hsa04630, FDR = 0.047), and extracellular matrix–receptor interaction (hsa04512, FDR = 0.05456). Cell adhesion molecules (hsa04514, FDR = 0.0055) and cytokines and inflammatory response (WP530, FDR = 0.137) were downregulated in the same cohorts 1 month after vaccination. However, delayed host immune response at Day 56 showed downregulation of NLR signaling pathway (hsa04621, FDR = 0.0029), PI3K‐Akt signaling pathway (hsa04151, FDR = 0.0057), and signaling by interleukins (R‐HSA‐449147, FDR = 0.0048), whereas inhibited IFN signaling (R‐HSA‐913531, FDR = 0.0014) and IFN‐α/β signaling (R‐HSA‐909733, FDR = 0.0015) were inactivated in patients at Day 84 postvaccination.

PPI network

On comparing combined postvaccination samples (Days 14, 28, 56, and 84) versus Day 0, 76 DEGs (41 up and 35 down) were defined. A PPI network was constructed from the combined DEGs. The clusters with densely connected nodes in the PPI network were detected (Figure 3G). KEGG enrichment analysis revealed 121 significant downregulated pathways. Eleven pathways were extracted in a separate subnetwork module; namely, influenza A, IL‐17 signaling pathway, TNF signaling pathway, chemokine signaling pathway, NLR signaling pathway, PI3K‐Akt signaling pathway, NF‐κB signaling pathway, JAK‐STAT signaling pathway, MAPK signaling pathway, Gap junction, and Leukocyte transendothelial migration (Figure 3H,I). The number of nodes was 2076, whereas the corresponding edges counted for 6123. To identify hub genes involved in the host immune response following BCG vaccination, a further GSEA of the cluster was performed. Downregulated genes with high strength of enrichment as NFKB1 (degree of connectivity with other genes = 91), RELA (degree = 95), IKBKB (degree = 60) genes were enriched in eight (out of 11) pathways, followed by PIK3CA (degree = 119), PIK3CD (degree = 59), PRKCB (degree = 26), RAF1 (degree = 69), AKT1 (degree = 221), and AKT2 (degree = 48) were involved in six (out of 11) pathways.

Comparison between transcriptomic signature in SARS‐CoV‐2 infection and BCG vaccination

Intersecting common pathways

A total of 52 enriched KEGG pathways were upregulated in COVID‐19 infection. On cross‐reference of these upregulated pathways with downregulated KEGG pathways in the BCG vaccination experiment, 45 common pathways were intersected accounting for 86.5% of SARS‐CoV‐2 upregulated pathways (Figure 4A–C). These pathways were categorized into the following groups: (1) Cellular processes, (2) organismal systems, (3) environmental information, and (4) human diseases. The top significant ones in the SARS‐CoV‐2 experiment were the IL‐17 signaling pathway (FDR = 1.39E−15, Hits = 14/92; overlapping genes in the pathway compared to the total gene set in the pathway), the TNF signaling pathway (FDR = 5.34E−15, Hits = 14/108 genes), and the NLR signaling pathway (FDR = 4.07E−10, Hits = 12/166 genes). Other immune‐related pathways were cytokine–cytokine receptor interaction (FDR = 2.08E−07, Hits = 11/263 genes), NF‐κB signaling pathway (FDR = 1.33E−07, Hits = 8/93 genes), and JAK‐STAT signaling pathway (FDR = 5.38E−06, Hits = 8/160 genes). In the same 45 pathways, viral infectious agents included influenza A (FDR = 1.10E−07, Hits = 10/168 genes), measles (FDR = 1.33E−07, Hits = 9/133 genes), and herpes simplex virus 1 infection (FDR = 1.33E−07, Hits = 10/181 genes), whereas bacterial infections included Legionellosis (FDR = 0.0014, Hits = 7/54 genes) and TB (FDR = 0.0039, Hits = 5/172 genes). Overlapping genes within each pathway that were significantly upregulated in SARS‐CoV‐2 infection and downregulated following BCG vaccination at the four different time points (14, 28, 56, and 84 days) are demonstrated in Table S13. Although different genes were modified at each time point, the pathways remained constantly downregulated post‐BCG vaccination.

Figure 4

Common pathways between SARS‐CoV‐2 and BCG vaccination experiments. (A) The intersection between KEGG pathways of SARS‐CoV‐2 infected cells and BCG vaccination. (B) The expression level of some hub genes following BCG vaccination showing downregulation postvaccination as an example of the reversed direction of expression. (C) Displays 45 upregulated KEGG signaling pathways in the SARS‐CoV‐2 cell line which are downregulated following BCG vaccination. Bars represent the observed gene count for each pathway. The degree of color represents the degree of significance (the more intensity, the higher significance). KEGG signaling pathways are categorized according to the functional hierarchical classification system. (D) The expression intensity of genes in the tuberculosis KEGG pathway. Colored by the log fold change of DEGs in (1) SARS‐CoV‐2 infection compared to mock‐treated cells and (2) following BCG vaccination. BCG, Bacille Calmette‐Guérin, KEGG, Kyoto Encyclopedia of Genes and Genomes

Comparative enrichment analysis in the TB pathway

On comparing 116 DEGs in SARS‐CoV‐2 infected cells which were over‐represented in the 179 gene set of TB pathways with their expression pattern after BCG vaccination, the top upregulated inflammatory markers (IL6/IL12, TLR1/2/9, TNF, C3, IL1A/B, NFKB1, TGFB2, AKT1, and HLA‐DPA1/B1), and transcription factors (NFKB1 and STAT1) showed reversed direction post‐BCG vaccination. Additionally, the downregulated DEGs in SARS‐CoV‐2, including several nonreceptor serine/threonine‐protein kinase (AKT3, CAMK2D, and MAPK10), a monocyte differentiation antigen (CD14), a membrane trafficking regulatory protein (LAMP2), growth factor (TGFB3), an apoptotic regulator (BCL2), and a proinflammatory cytokine (IL18) switched their direction following BCG vaccine administration (Figure 4D). Similar paradox directions of expression level in some DEGs for SARS‐CoV‐2 cells and BCG vaccine samples were also noted in Figure S4–S10.

DISCUSSION

Overview

To gain insight into the engendering mechanism of the SARS‐CoV‐2 infection, the gene expression profiles of the virus were systematically analyzed through bioinformatics techniques. In this study, a total of 161 DEGs, including 113 upregulated and 48 downregulated, were screened. The biological functions of these DEGs were explored based on GO function and pathway enrichment data. Further analysis indicated significant upregulation in 45 distinct KEGG pathways in SARS‐CoV‐2 infection overlapped with pathways downregulated following the BCG vaccination. In addition to the pathway analysis, this study describes potential kinases, transcription factors, miRNAs, and lncRNAs differentially expressed in Sars‐CoV‐2 infection (Figure 2C). Several candidate DEGs, including IL‐6, CCL20, CSF2, ICAM, and CXCL1/2 were highlighted for their key roles in respiratory viral infection, TB, and overall immune function. To elucidate the relationship between COVID‐19 infection and BCG‐vaccination, the common pathways were categorized into the following groups: Inflammatory and immunoregulatory, signaling, and infectious pathways.

Common inflammatory and immunoregulatory pathways

Of SARS‐CoV‐2 infected cells, analysis implicates several upregulated genes enriched in inflammatory pathways, most significantly in the IL‐17 and NLR pathways. Conversely, these same pathways were downregulated in the BCG‐vaccinated group, suggesting BCG vaccination may compensate for aspects of pathway dysregulation in SARS‐CoV‐2 infection. In IL‐17 signaling, the family of IL‐17 mediates protective innate immunity against external pathogens and plays a central role in the self‐clearance of intracellular pathogens. , Additionally, T helper cells (Th17), which themselves produce Il‐17, are known to be key players in the pathogenesis of chronic inflammatory diseases and autoimmune tissue destruction. Elevated Th17 responses and IL‐17 pathways are seen in COVID‐19 patients and have been linked to “cytokine storms,” a surge of proinflammatory molecules that are associated with the acute respiratory distress syndrome (ARDS) typically seen in these patients. An antagonistic effect is noted in BCG vaccination, which has been shown to decrease Th17 cell maturation, suppress Th17 response, and prevent the production of IL‐17 in the Rhesus macaque model, thereby conferring nonspecific immunity. , Furthermore, the IL‐17 gene family eliminates self‐reactive T cells through the production of granulocyte colony‐stimulating factor (G‐CSF) and several chemokines upregulated in Sars‐CoV‐2 infected cells, including CXCL1, CXCL2, CCL20, and IL‐17 , (Figure 2D).

The other major upregulated inflammatory pathway, the NLR‐signaling pathway, detects various pathogens and stimulates innate immune responses against them, driving the subsequent activation of cytokine, NF‐κB, and MAPK pathways. Upregulation of NLR signaling without appropriate negative feedback regulation can contribute to pathological tissue damage. Innate and adaptive immune response to an invading pathogen relies on the ability of the body to recognize foreign elements and is stimulated by inflammasomes, which are intracellular multiprotein complexes such as NLR.  NLRC5 (short for NLR family CARD domain containing 5) regulates major histocompatibility complex class I expression during a viral infection. In accordance with this, a study reported knockdown of NLRC5 resulted in decreased levels of CD8+ cells in influenza A virus. Like IL‐17 signaling, the upregulation of NLR signaling may further be contributing to the cytokine storms evident in COVID‐19 patients. These same pathways are markedly downregulated post BCG vaccination, suggesting BCG vaccination could potentially attenuate these pathways.

Common signal transduction pathways

Several signal transduction pathways were upregulated in Sars‐CoV‐2 and inversely downregulated in the BCG‐vaccinated group with TNF, NF‐κB, MAPK, and JAK/STAT signaling pathways conferring the most statistical significance. TNF signaling upregulation is heavily involved in the complex regulation of immune cells. Following trimerization with either TNFR1 (expressed nearly ubiquitously) or TNFR2 (expressed mostly in immune cells), the TNF signaling pathway can activate NF‐κB and MAPK signaling pathways downstream of it. Upon viral infection, NF‐κB proteins regulate the transcription of many genes, including antimicrobial peptides, cytokines, chemokines, stress‐response proteins, and antiapoptotic proteins. In a study of Middle Eastern Respiratory Syndrome, a viral infection caused by a member of the coronavirus family, researchers found the virus downregulated antiviral cytokines (TNF‐α) as it induced proinflammatory cytokines (IL‐1β, IL‐6, and IL‐8) during initial infection. In a study of SARS, yet another member of the coronavirus family, TNF‐α was similarly downregulated and induced NF‐κB activation. In the SARS pathogen, ACE2 was determined to be the functional receptor of the virus through the regulation of TNF‐α, which ultimately permitted viral entry and promoted respiratory pathogenesis. Altogether, the literature on other members of the coronavirus family indicates disruption of the TNF‐signaling pathway, which is in accordance with the marked upregulation of TNF pathway signaling noted in our findings. In BCG vaccinated models in vitro, exposure of human peripheral blood mononuclear cells (PBMCs) to BCG treatment boosted IL‐6 and TNF‐α expression in response to lipopolysaccharide stimulation. Similarly, BCG‐immunized adults produced high TNF‐α and IL‐1β expression 3‐month postvaccination. In addition to sustained cytokine response at 3 months, the production of heterologous Th1 and Th17 remained elevated 1‐year post BCG vaccination, concluding the vaccine created long‐term immune responses to pathogens besides M. tb alone. Again, we note the same TNF pathway upregulated in the SARS‐CoV‐2 virus is downregulated in BCG vaccination. Downstream of TNF, the NF‐κB signaling pathway similarly affects a broad range of biological processes, including adaptive immune, inflammatory, and stress responses. Upon viral infection, NF‐κB proteins regulate the transcription of many genes, including antimicrobial peptides, cytokines, chemokines, stress‐response proteins, and antiapoptotic proteins. Among the SARS‐CoV‐2 upregulated genes in the NF‐κB pathway were CSF‐2 and IL‐6 are of particular note given their role as innate immune system mediators. , The upregulation of CSF‐2, produced by endothelial and immune cells, is similarly upregulated in other respiratory diseases, including mycoplasma pneumoniae and M. tb, by promoting neutrophil and macrophage inflammatory response. The upregulation of IL‐6 has been noted in other members of the coronavirus family, including the pathogenic agent responsible for SARS. In murine models, SARS‐coronavirus spike protein, which is evolutionarily conserved in Sars‐CoV‐2, induced the upregulation of IL‐6 through the NF‐κB pathway. Higher levels of proinflammatory cytokines, including IL‐6, TNF, CSF‐2, IL‐1b, and IL‐8 are noted in COVID‐19 patients with increasing levels as a predictor of the severity of pneumonia.

Additionally, MAPK signaling was upregulated in the Sars‐CoV‐2 infected group but inversely downregulated in the BCG‐vaccinated group. MAPK regulates cellular processes of proliferation, stress responses, as well as immune defense. , ,  The MAPK signaling cascade is activated in response to external stress signals of the three MAP kinases (ERK, JNK, and p38 isoforms) with subsequent stimulation of proinflammatory cytokines by JNK and p38 signaling. , , In the nucleus, activated JNK and p38 promote multiple effector proteins, including NF‐κB, c‐Jun, STAT1. These effector proteins have been implicated in viral infections such as influenza A and HSV‐1, and their dysregulation by pathogens is associated with an impaired antiviral response by the host. , The downregulation of the MAPK pathway in the BCG‐vaccinated group (Figure 3I) may be able to partially antagonize the upregulation induced by the SARS‐CoV‐2 virus. Additionally, our results indicated the top three upregulated kinases in the SARS‐CoV‐2 group were MAKP14, MAPK3, and MAPK1, further corroborating the central role MAPK plays in viral pathogenesis of SARS‐CoV‐2 infection (Figure S6).

Lastly, JAK‐STAT signaling, which is a major pathway involved in cytokine signaling and subsequent inflammation, was upregulated in SARS‐CoV‐2 infection as well. IFNs induce the JAK/STAT signaling pathway, which in turn activates transcription of widespread immune protection. For this reason, several viruses have evolved to target the JAK/STAT pathway themselves. Several studies have recently advocated for the potential benefit of commercially available Jak 1 and 2 inhibitors as an anti‐inflammatory treatment in COVID‐19 cases ; however, there is concern from the scientific community that JAK inhibitors may promote the evolution of SARS‐CoV‐2 virus as this has been reported in herpes viruses. BCG vaccination may be key in attenuating the JAK/STAT pathway without entirely impairing IFN‐mediated response. In vitro experiments show that BCG vaccination induces two members of the suppressor of cytokine signaling (SOCS) family, SOCS1/3, eliciting a negative feedback regulator of the JAK/STAT signaling cascade via IFN‐γ regulation.

Overall, although TNF, NF‐κB, MAPK, and JAK/STAT signaling were upregulated in SARS‐CoV‐2 infection, they were markedly downregulated in BCG, suggesting BCG vaccination may be targeting these same pathways. The exact mechanism remains elusive, but identification of these pathways and their targets is key to mitigating pathogenesis and merits further investigation.

Infectious and human disease‐signaling pathways

BCG vaccination downregulated several pathways related to viral and bacterial inflammatory disease pathways as well, including influenza A, measles, herpes simplex virus 1, legionellosis, and TB, all of which were upregulated in SARS‐CoV‐2 infection. Prior studies have demonstrated nonspecific effects of BCG against viral infection and bacterial infections. , In one such study, mice inoculated with BCG displayed overall increased resistance to encephalomyocarditis, murine hepatitis, type 1 and 2 herpes simplex, foot‐and‐mouth disease, and A0 and A2 influenza viruses. Intercellular adhesion molecules (ICAM) included in the upregulated DEGs (Table S1), are receptors exploited by viruses for entry into host cells, intercellular signaling, and continued survival. In the human rhinovirus and influenza virus, ICAM‐1 is involved in viral protein uncoating, delivering the viral RNA genome into the cytoplasm of host cells across the lipid bilayer. , Other studies corroborate this, showing that most viruses interact with and induce ICAM‐1 expression, including HIV, human parainfluenza virus, and rhinovirus infection. , , The presence of cytokines such as IFN‐γ, TGF‐β, and TNF‐α induces the expression of ICAM on human TB infected macrophages, further activating T‐Cells. Intercellular adhesion molecules (or CD54), included in the upregulated DEGs, are receptors utilized for pathogen entry into host cells that promote intercellular signaling and survival of the virus in the host cell. In human rhinovirus infection and influenza virus, ICAM‐1 is involved in viral protein uncoating, delivering the viral RNA genome into the cytoplasm of host cells across a lipid bilayer. , Other studies corroborate this showing that most viruses interact with and induce ICAM‐1 expression, including HIV, human parainfluenza virus, and rhinovirus infection. , , The presence of cytokines such as IFN‐γ, TGF‐β, and TNF‐α induces expression of ICAM on macrophages for activation of T cells, this was observed in human TB infected macrophages. BCG immunotherapy of bladder tumor cells showed mycobacteria infection‐induced expression of ICAM‐1 molecules, thereby eliciting immune response through improved antigen presentation to T lymphocytes. BCG is thought to inhibit inflammasome activation via zinc metalloproteases to improve immunogenicity , and remains the current standard therapy for non‐muscle‐invasive bladder cancer.

The role of VDR and other enriched transcription factors

Compared to mock‐treated cells, SARS‐CoV‐2 infected cells exhibited significant enrichment of VDR transcription factor, one of the major downstream signaling transduction nuclear receptor in TB KEGG pathway, and was significantly associated with the upregulated DEGs, namely, NFKB2, NFKBIA, and IRAK2, which are related to NF‐κB and Toll‐like receptor signaling pathways. Normally, VDR binds to calcitriol, the active form of vitamin D. The VDR‐calcitriol complex then interacts with retinoid‐X‐receptor to form a heterodimer capable of regulating transcriptional activity, leading to several pleiotropic transcriptional effects. Similar to its upregulation in SARS‐CoV‐2, VDR gene expression is continuously upregulated during HIV and influenza virus infection. Previously reported, VDR signaling may repress cytokine gene expression in activated T‐cells, consequently reducing inflammatory response. , Part of the mechanism of action of the BCG vaccine against TB is through increased production of vitamin D, which binds with its nuclear receptor VDR, resulting in the generation of antimicrobial peptides (cathelicidin and β‐defensin) and death of intracellular M. tb. In short, vitamin D plays an important role in modulating the innate and adaptive immune response alongside its classically characterized role in bone health. Long before the scientific basis was understood, vitamin D was unknowingly and empirically being used to help TB patients who were sent to sanatoriums for sun‐light exposure and prescribed cod‐liver‐remedies rich in vitamin D. Since then, vitamin D has been well characterized for its immunomodulatory properties, including the regulation of IFN‐γ, a key activator of macrophage response, and promoting the production of regulatory T cells. In addition to its role in M. tb, cathelicidin is a highly conserved protein that has been shown to direct antiviral activity in several respiratory viral infections including respiratory syncytial virus, human rhinovirus (HRV), and influenza A. , , , , ,  In vitro, vitamin D has been shown to counteract the M. tb‐induced downregulation of cathelicidin by actually recovering cathelicidin levels and promoting Th1 cell differentiation. Following BCG vaccination, vitamin D levels remain elevated, indicating vitamin D is being upregulated long after the initial inoculation. It has been postulated that increased calcitriol recruits dendritic cells from the site of inoculation to the lymph nodes, accounting for the sustained response. This serves as a possible mechanism of protective effects of BCG vaccination in relation to consistently elevated plasma vitamin D concentration present for as long as 9 months as reported by the authors. Vitamin D deficiency has been shown to contribute to the pathogenesis of ARDS. In children, the mutation in the VDR gene single‐nucleotide polymorphism reported in children was associated with a viral infection. Given these findings, the transcriptional upregulation of VDR in SARS‐CoV‐2 may be contributing to the dysregulation of immunological response by upregulating cytokine response, interrupting microbial peptides, and promoting immune invasion in the lungs; therefore, the well‐characterized role of VDR in BCG vaccination could serve as a potential biomarker against COVID‐19 disease.

In addition to VDR, our bioinformatic analysis showed the deregulation of multiple transcription factors (TFs) in SARS‐CoV‐2 infection. Specifically, TAL2, DACH2, and AFF2 have upregulated; however, the functional role in the viral immune response is unclear. Conversely, transcription factors HEY2, KLF2, EGR2, and ZBTB16 were under‐expressed and were previously reported to have an eminent role in the fine‐tuning immune response. Among these, HEY2 was the most downregulated TF gene. HEY2 is a known transcription repressor, modulating cardiovascular development, neurogenesis, and oncogenesis. A recent study demonstrated the role of HEY2 in inflammation, namely, chronic periodontitis, via negative regulation of IL‐6, IL‐1β, and TNF‐α expressions. Similarly, KLF2 is a key player in T cell differentiation, trafficking, quiescence, and survival, in addition to the regulation of endothelial function. Mice with KLF2 deficiency experienced an activated T cell phenotype with severely reduced T cells in the periphery and increased susceptibility to HIV‐1 infection. , Next, the EGR2 transcription factor can control adaptive innate immunity and lead to functional impairment of T cells. After infection with influenza virus, EGR2 knockout mice exhibited prolonged viral shedding, impaired CD4+ T‐cell response, and infiltration of memory precursor type CD8+ T cells into the lung with decreased IFN‐γ, TNF‐α, and granzyme B. Lastly, the ZBTB16 transcription factor can mediate the ubiquitination and subsequent proteasomal degradation of target proteins. Transcription factor enrichment analysis unravels multiple transcription factors modulating the deregulated genes in SARS‐CoV‐2 infected cells. RELA, a well‐known antiviral transcription factor, was the top enriched transcription factor for DEGs following SARS‐CoV‐2 infection. It is previously reported to promote the growth of cytopathic RNA viruses by extending the lifespan of infected cells and serve as the replicative niche of intracellular pathogens.  TP63 was also shown to be a cellular regulator of the human papillomavirus life cycle.

MicroRNAs, snoRNAs, and pseudogenes

In the present study, we observed some deregulated miRNAs and highlighted the relationship between DEGs and miRNA networking which may account for the transcriptomic changes evident in SARS‐CoV‐2 expression. Of particular note, mir‐26a‐5p, miR‐26b‐5p, and miR‐124‐3p, which have all been well characterized for their involvement in viral and bacterial inflammatory pathways, particularly influenza A, RSV, and M. tb, were noted here for their putative respiratory pathogenesis. In M. tb infection, miR‐26a‐5p can modulate macrophage IFN‐γ responsiveness. In another study, this same microRNA was shown to be downregulated in influenza A virus in humans, contributing to viral pathogenesis. The other family member, miR‐26b‐5p, can inhibit viral replication in the vesicular stomatitis virus and Sendai virus by inducing type‐1 IFN expression. RSV acts to upregulate miR‐26b, which in turn inhibits Toll‐like receptor TLR4, a key sentinel in adaptive immune response. A similar effect on TLRs is seen in the pathogenesis of M. tb in alveolar macrophages wherein microRNA‐124 negatively regulates TLR signaling. SARS‐CoV‐2 infected cells exhibited overexpression of miR‐23a, which was previously shown to promote the replication of human herpes simplex virus type 1 via downregulating interferon regulatory factor 1 (IRF1), an innate antiviral molecule. Altogether, the highlighted cluster of microRNAs is likely regulating viral pathways relevant to SARS‐CoV‐2 infection. Posttranscriptional regulation of some genes altered in SARS‐CoV‐2 infected cells are, therefore, worthy of experimental validation.

snoRNAs and pseudogenes were found to be deregulated in SARS‐CoV‐2 cells. SnoRNAs can serve as a source of short regulatory RNA species involved in the control of processing and translation of various messenger RNAs. Despite their exact role in COVID‐19 and TB being unclear and complicated, these ncRNAs were found to be upregulated in virus‐infected human cells. They can guide chemical modifications of structural RNAs and have a role in cell–cell communication. SnoRNAs were reported to have a paradox effect; they can act as mediators of host antiviral response, and on the contrary, are utilized by viruses to evade innate immunity and complete their life cycle. Furthermore, different groups of pseudogenes were activated in human cells upon triggering inflammation with known cellular protein, including virus particles and chunks of bacterial cell walls, highlighting the unique functions of pseudogenes in host immune response.

LIMITATIONS

To the best of our knowledge, our study is the first to inquire into the regulatory mechanisms underlying the host immune response against SARS‐CoV‐2 infection and to provide a putative protective role of BCG vaccination. Altogether, several DEGs with a high connectivity degree in the PPI network were identified in this study, which showed the strong inverse correlation in SARS‐CoV‐2 infection and BCG vaccination with upregulation in the former and downregulation in the latter. This inverse relationship suggests BCG vaccination may mediate key pathways in SARS‐CoV‐2 infection through NSEs. Nevertheless, it is essential to recognize some limitations of administering BCG vaccines, including, the heterogeneity of BCG response due to pharmacogenomic variation of patients, contraindications of BCG‐vaccination for immunosuppressed patients, and BCG response variations in age groups. Although the exact mechanisms with which BCG elicits protective effects are not entirely understood, its efficacy in infectious, oncogenic, and inflammatory diseases are well documented and may confer similar effects in COVID‐19. More research is warranted to determine whether the pathways upregulated genetically in SARS‐CoV‐2 infected cells function the same in vivo.

CONCLUSION

In summary, the bioinformatics approach identified key DEGs involved in the pathogenesis of COVID‐19. Significant nodes in the PPI network, including high degrees in IL‐17 signaling, TNF signaling pathway, NLR, and NF‐κB give novel insights into the inflammatory and signal transduction pathways mediating SARS‐CoV‐2 viral infection. The 45 common pathways upregulated in SARS‐CoV‐19 and downregulated in BCG vaccination suggest BCG vaccination may help to mitigate this pathway dysregulation. Although further experimental validation is warranted to verify our discoveries, the results give insights into potential biomarkers and targeted therapy to address the ongoing COVID‐19 pandemic. We maintain BCG vaccination may be a safe and cost‐effective alternative in incurring partial protection against the COVID‐19 pandemic until more targeted measures can be produced and implemented.

CONFLICT OF INTERESTS

The authors declare that there are no conflict of interests.

AUTHOR CONTRIBUTIONS

Eman A. Toraih and Emad Kandil designed the study. Eman A. Toraih and Mohammad H. Hussein analyzed the data. Eman A. Toraih and Emad Kandil interpreted the data. Eman A. Toraih, Jessica A. Sedhom, and Titilope M. Dokunmu searched the literature and wrote the manuscript. All authors reviewed and approved the final version of the manuscript.

Supporting information.

Click here for additional data file.

Supporting information.

Click here for additional data file.