Identifying molecular insight of synergistic complexities for SARS-CoV-2 infection with pre-existing type 2 diabetes.

The ongoing COVID-19 outbreak, caused by SARS-CoV-2, has posed a massive threat to global public health, especially to people with underlying health conditions. Type 2 diabetes (T2D) is lethal comorbidity of COVID-19. However, its pathogenetic link remains unclear. This research aims to determine the genetic factors and processes contributing to the synergistic severity of SARS-CoV-2 infection among T2D patients through bioinformatics approaches. We analyzed two sets of transcriptomic data of SARS-CoV-2 infection obtained from lung epithelium cells and PBMCs, and two sets of T2D data from pancreatic islet cells and PBMCs to identify the associated differentially expressed genes (DEGs) followed by their functional enrichment analyses in terms of protein-protein interaction (PPI) to detect hub-proteins and associated comorbidities, transcription factors (TFs), microRNAs (miRNAs) as well as the potential drug candidates. In PPI analysis, four potential hub-proteins (i.e., BIRC3, C3, MME, and IL1B) were identified among 25 DEGs shared between the disease pair. Enrichment analyses using the mutually overlapped DEGs revealed the most prevalent GO and cell signalling pathways, including TNF signalling, cytokine-cytokine receptor interaction, and IL-17 signalling, which are related to cytokine activities. Furthermore, as significant TFs, we identified IRF1, KLF11, FOSL1, and CREB3L1 while miRNAs including miR-1-3p, 34a-5p, 16-5p, 155-5p, 20a-5p, and let-7b-5p were found to be noteworthy. The findings illustrated the significant association between COVID-19 and T2D at the molecular level. These genetic determinants can further be explored for their specific roles in disease progression and therapeutic intervention, while significant pathways can also be studied as molecular checkpoints. Finally, the identified drug candidates may be evaluated for their potency to minimize the severity of COVID-19 patients with pre-existing T2D.

Introduction

After the outbreak of SARS-CoV [1] and MERS-CoV [2] in 2002 and 2012, respectively, a novel coronavirus SARS-CoV-2 has emerged in Wuhan city of Hubei province of China at the end of 2019 with its massive infectious threat [3,4]. This virus causes respiratory tract infection with clinical syndromes including fever, cough, sore throat, pneumonia, and in some severe cases, acute respiratory distress syndrome (ARDS), sepsis and septic shock, multiorgan failure, such as acute kidney injury and cardiac arrest [5,6]. The virus SARS-CoV-2 has been identified as the causative pathogen of the outbreak of coronavirus disease (COVID-19) by deep sequencing and etiological investigations [7]. The World Health Organization (WHO) has already declared COVID-19 as a pandemic. As of June 29, 2021, 181,007,816 cases have been confirmed as SARS-CoV-2 infected in more than 218 countries and regions with 3,927,222 deaths according to WHO [8].

Comorbid diseases, such as diabetes, hypertension, and heart diseases have been reported as risk factors for COVID-19 patients [6,9]. In Ref. [6], it was reported that, out of 41, 20 SARS-CoV-2 infected patients had comorbidities including diabetes (20%), hypertension (15%), and cardiovascular disease (15%). In Ref. [9], 191 COVID-19 patients’ records from two hospitals of Wuhan, China were investigated. Among them, 91 (48%) patients had comorbidity, of which 58 (30%) patients had hypertension, 36 (19%) had diabetes, and 15 (8%) were suffering from coronary heart disease. Moreover, several recent studies reported that type II diabetes (T2D) increases the risk of infection, severity, and mortality of COVID-19 patients [10,11]. However, there is a lack of sufficient information regarding genetic interactions of T2D on COVID-19.

Viruses are obligate intracellular pathogens and they cannot replicate without host cellular factors during infection [12]. Consequently, virus-host protein-protein interaction (PPI), might be an effective way toward elucidating the mechanism of viral infection [13]. Recently, integrative network-based approaches showed better effectiveness to identify disease-associated biomarkers and therapeutic targets [14] of various diseases, including respiratory system diseases [15], neurodegenerative diseases [[16], [17], [18]], cardiomyopathy [19], and viral infections, such as SARS-CoV, HIV [20] and Zika [21]. The network-based strategies can potentially be used to identify effective repurposable drugs [[22], [23], [24]] as well as the combination of drugs [25] for various human diseases as well. Recently, genetic investigations into the gene expression data offered a better understanding of the molecular mechanism of the SARS-CoV-2 infection and its genetic association with various complications [[26], [27], [28], [29], [30]].

In this study, we investigated the genetic interaction of SARS-CoV-2 with T2D using this network-based strategy that incorporates studies regarding gene expression profiling, protein-protein interaction, gene ontologies, molecular pathways, and regulatory analysis. The findings may lead to determining the significant therapeutic target to fight against the ongoing pandemic due to SARS-CoV-2 infection. The analytical approach adopted in this research work is illustrated in Fig. 1 .

Fig. 1

Schematic view of the systematic pipeline used in this study. (A) At first, transcriptomic analyses of two RNA-Seq data, one NanoString and one microarray data and cross-comparison identified 15 common DEGs between SARS-CoV-2 infection and T2D for lung epithelium and pancreatic islet cells respectively, and 11 common DEGs for peripheral blood mononuclear cells (PBMCs). (B) Biological functions of these 26 (25 distinct) DEGs were assessed by PPI analysis and functional enrichment analysis using GO and cell signalling pathway databases. (C) Therapeutic targets were identified by obtaining hub genes and putative drug candidates. (D) Regulatory elements and possible comorbidities were determined. (E) All the gained results were validated through an extensive literature review.

Materials and methods

Data

To reveal the genetic association of SARS-CoV-2 with T2D, we collected the RNA-Seq and gene expression microarray data from National Center for Biotechnology Information Gene Expression Omnibus (NCBI-GEO) and EBI array express. In this study, we analyzed the datasets with accession numbers GSE147507 and E-MTAB-8871 for SARS-CoV-2 while the datasets with accession numbers GSE9006 and E-MTAB-5060 for T2D. The GSE147507 is an RNA-Seq data obtained by transcriptional profiling of SARS-CoV-2 infected and mock-treated cells from humans and ferrets using Illumina NextSeq 500 platform [31]. Herein, the human lung epithelium (NHBE) cell replicates from 3 SARS-CoV-2 and 3 mock-treated samples were collected. The E-MTAB-8871 was obtained from human peripheral blood mononuclear cells (PBMCs) of 10 healthy individuals and 23 samples from 3 SARS-CoV-2 infected patients using NanoString profiling [32]. The T2D dataset (GSE9006) was obtained by gene expression microarray using Affymetrix HG-U133 array of PBMCs collected from 24 healthy individuals and 12 T2D patients [33]. The E-MTAB-5060 was produced by RNA-sequencing of human pancreatic islet cells from 4 T2D patients and 3 healthy individuals using Illumina TruSeq protocol [34]. Since the complications and severity of COVID-19 usually vary drastically, we compiled two sets of experimental data to ensure the validity of our results. Firstly, we analyzed the differentially expressed genes (DEGs) for SARS-CoV-2 infection and T2D with data obtained from lung epithelial and pancreatic islet cells, respectively. Secondly, we compared the DEGs for the same disease-pair with data collected from human PBMCs.

Identification of significant dysregulated genes

The RNA-Seq is a next generation sequencing technology to measure gene expression with a high level of accuracy and mitigates many limitations of microarrays [35]. Using this high-throughput sequencing technology and global transcriptome analyses, we compared the gene expression profiles of SARS-CoV-2 and T2D. To identify DEGs associated with the respective diseases, we used an R Bioconductor package DESeq2 [36] which is explicitly designed for RNA-Seq data. This package uses the Student unpaired -test for the identification of DEGs.

As the GSE9006 dataset is a microarray data set, we analyzed it with the R Bioconductor package Limma [37] to obtain the dysregulated genes. Since various errors can be introduced in preparing and analyzing microarray data of different platforms and experimental systems, the gene expression data in each sample (case or control) were normalized using the -score transform which is defined aswhere and are the mean and standard deviation of the expression value of -th gene estimated over all the samples, respectively. After obtaining the DEGs for each disease condition, we selected the significant genes by setting the threshold level of the absolute value of base-2 log Fold Change 1 and adjusted -value (FDR) . The shared DEGs by the SARS-CoV-2 and T2D were then obtained through a cross-comparative analysis.

Profiling associated comorbidities

To assess the possible health implications linked to COVID-19, we anticipated the comorbidities associated with the shared DEGs via Enrichr web-platform [38] using DisGeNET database [39]. This database is a publicly available repository of genes that are associated with human diseases and the latest version (v7.0) includes information regarding the association of 21,671 genes with 30,170 disorders (https://disgenet.org). In this analysis, we only considered diseases or complications with gene enrichment and adjusted -value. Likewise, we ran the shared DEGs using the same database via Metascape server [40], where statistical significance was set to adjusted p-value. For visualization of the gene-disease association (GDA), a bipartite network was created and designed with Cytoscape v3.8 software [41]. It is an open source software to produce high resolution interactomes and complex networks.

Functional enrichment analysis

Functional enrichment analysis, also referred to as gene set enrichment analysis (GSEA), identifies a group of significantly enriched genes engaging in the course of biological action or residing at the central chromosomal location by integrating all prior knowledge [42]. Such a group of genes can be associated with a biological function by means of gene ontology (GO) that has 3 separate classes including biological process (BP), molecular function (MF), and cellular component (CC) [43]. Similarly, pathway analysis also plays an important role to reveal the molecular or biological function that underlies the development of complex diseases and it also helps us to understand how to intervene therapeutically in disease processes. Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway database constitutes metabolic pathways representing the network of molecules interconnected for a general function [44]. Besides KEGG, Wikipathways [45], and Reactome [46] pathway databases are also widely used as annotation sources for analyzing signalling pathways. We performed the GSEA using the online platform Enrichr by integrating the pathway and GO term annotations for all the shared DEGs between the SARS-CoV-2 and T2D to have a better understanding of the metabolic pathways active in SARS-CoV-2. Enrichr integrates available gene-set libraries to provide enrichment outcomes for a gene list of interest. For statistical significance, the adjusted -value was considered as less than 0.05 for the significance assessment of the obtained enrichment results.

Protein-protein interaction analysis

Proteins usually interact with each other to form molecular machines that perform cellular functions. Such protein-protein interactions (PPIs) can be mapped to discern the functional and structural knowledge of cellular protein networks [47]. We incorporated the web-based visualization software NetworkAnalyst v3.0 to construct the PPI network of the common DEGs. This web-tool interprets the gene expression data from PPI network perspective [48]. Herein, we employed the STRING database with a confidence score of 600 for network construction. This database comprises published and estimated interactions that include both functional and physical contacts [49]. In this case, we only considered an experimentally validated interaction dataset. To determine potential hubs within the PPI network, we then applied three different methods, i.e., degree, maximal clique centrality (MCC), and betweenness using cytoHubba plugin [50] in Cytoscape v3.8. Next, we compared the results and identified the common nodes as the most potential hubs. Finally, the obtained networks were customized in Cytoscape v3.8.

Probing the regulatory networks

Transcription factors (TFs) control the gene expression level by switching it on and off at the transcription level. Thus, TFs and their interactions with genes are essential for their functional implications [51]. We utilized the NetworkAnalyst platform to investigate the interactions of TFs with the shared DEGs. For this purpose, we selected the Encyclopedia of DNA Elements (ENCODE) database in NetworkAnalyst that uses peak intensity signal 500 and the predicted regulatory potential score 1. The ENCODE project provides genomic and transcriptomic annotations to characterize the human genomic elements [52]. Likewise, microRNAs (miRNAs) are non-coding RNAs that affect gene expression via post-transcriptional regulation [51]. To generate the gene-miRNA interactome, we provided the shared DEGs into the same platform utilizing the TarBase 8.0 database, which provides an experimentally validated collection of curated miRNA-gene interaction data as well as function-related data derived from various present-day empirical methods [53]. Both gene-miRNA and gene-TF networks were initially filtered with the degree centrality. Next, the obtained networks were analyzed with degree method via CytoHubba plugin to determine the key regulatory components and their interacting gene partners.

Identification of drug candidates

One of the primary objectives of this line of research focuses on pinpointing the potential drug molecules. We searched the Drug Signatures Database (DSigDB) based on the shared DEGs between SARS-CoV-2 and T2D using the Enrichr web-platform. The DSigDB is a collection of 22,527 gene sets related to the drug and small molecules considering the dysregulation in gene expression due to drug/compounds [54]. The significant drug molecules were filtered through manual curation based on adjusted -value.

Results

Shared differentially expressed genes

We first compared the gene expression levels of lung epithelium cells from SARS-CoV-2 patients (GSE147507) and pancreatic islet cells from T2D individuals (E-MTAB-5060) with the corresponding cells from healthy controls. Similarly, we also compared the expression profiling of PBMCs from SARS-CoV-2 (E-MTAB-8871) and T2D (GSE9006) patients with normal individuals. By filtering DEGs for each disease condition, we found 156 and 78 significant DEGs in the first experiment for SARS-CoV-2 and T2D, respectively. In the second experiment, we found 1289 and 73 significant DEGs in PBMCs for the same disease-pair. The four volcano plots in Fig. 2 presented the expression pattern of the genes in both experiments.

Fig. 2

Differential gene expression and common DEGs. Volcano plots depict the genes expression in A) SARS-CoV-2 infected lung epithelium cells, B) SARS-CoV-2 infected PBMCs, C) PBMCs of T2D patients, and D) T2D diseased pancreatic islet cells (red dots indicate significant DEGs), while bubble plot shows (E) the expression pattern of 25 common DEGs between SARS-CoV-2 and T2D in lung epithelium cells and pancreatic islet cells, respectively as well as in PBMCs.

The cross-comparative analysis revealed 15 common DEGs between SARS-CoV-2 and T2D for lung epithelium and pancreatic islet cells, respectively. On the other hand, for PBMCs, we found 11 common DEGs between the disease pair. By combining these two sets of DEGs, we found a total of 25 unique shared DEGs. Negative 10-base logarithmic transformed FDR adjusted -value of these 25 common DEGs are shown in Fig. 2-E. Notably, all the common genes were strikingly up-regulated for the lung epithelium and pancreatic islet cell in SARS-CoV-2 and T2D, respectively, whereas CXCL1 was irregularly expressed in both scenarios.

Significant gene ontologies and signalling pathways

After identifying the overlapping DEGs for the disease pair, we performed extensive GO and signalling pathway analyses accessing relevant curated databases. The adjusted -value filtration found a total of 38 GO terms (Fig. 4) and 46 molecular pathways (Fig. 5) to be significantly enriched. Curation found cytokine activity and cytokine-mediated signalling pathway to be predominant. Similarly, signalling pathway analyses exhibited TNF signalling pathway and cytokine-cytokine receptor interaction having noteworthy enrichment.

Fig. 4

Gene ontology analysis revealed significant GO terms associated with SARS-CoV-2 and T2D. The biological process, cellular component, and molecular function datasets were considered for this analysis.

Fig. 5

Significant signalling pathways associated with SARS-CoV-2 and T2D. The human KEGG (2019), reactome (2016), and human WikiPathways (2019) datasets were considered for this analysis.

Genetic relationship of COVID-19 with other diseases

By analyzing the GDA, we identified 19 highly significant (-value ) diseases associated with up to 18 shared DEGs. Fig. 3 shows the different cancers and diseases related to DEGs identified in COVID-19. Most of the comorbidities were different forms of cancers. The highest number of DEGs (18) were found to be associated with breast cancer and its metastasis while at least 10 DEGs were affiliated with ulcerative colitis, cervix carcinoma, bladder neoplasm, atherosclerosis, autoimmune diseases, and inflammatory bowel disease (IBD). Other diseases include rheumatoid arthritis (13), asthma (12), Alzheimer's disease (11), and arteriosclerosis (11). Additionally, the Metascape server predicted the top 20 diseases that are associated with 8–11 genes as shown in Fig. 3-B. These DEGs and diseases could provide the directionality of the COVID-19 progression and its long-term health impact.

Fig. 3

Gene-disease association network. In this figure, (A) the bipartite network includes circular nodes (blue) representing the shared DEGs and hexagonal nodes indicating COVID-19 (yellow) and different diseases (green), and (B) the bar graph indicating the top 20 diseases associated with DEGs in terms of expression and number of DEGs involved as predicted by the Metascape server.

Hub-proteins identified in the PPI network

We constructed the PPI network using all the distinct DEGs shared between the disease pair through NetworkAnalyst. The resulted PPI network comprises 171 nodes and 176 edges as shown in Fig. 6 -A. Each node in the network represents a protein and an edge indicates the functional interaction between two proteins. We employed three methods for topological analysis where each method identified the top 7 hub-nodes within the PPI network (Fig. 6B–D). Interestingly, four out of 7 hub-proteins were anticipated by all methods. As validated by multiple methods and exhibited a minimum of eight interconnections, we recognized those four hub-nodes as potential hub-proteins, i.e., BIRC3, C3, MME, and IL1B.

Fig. 6

The protein-protein interaction network and hub-proteins. This network was constructed with the DEGs shared by SARS-CoV-2 and T2D using STRING database (confidence cut off 600). The network depicting (A) a total of 171 proteins including 11 shared DEGs in which 7 hubs were indicated as predicted by the degree method. Additionally, three smaller networks are depicting hub-proteins anticipated by (B) degree, (C) betweenness, and (D) maximal clique centrality (MCC) methods. For all methods, the top seven hub-proteins are indicated by the color ranging from red (higher) to yellow (lower).

Regulatory networks revealed highly active TFs and miRNAs

The TF-gene interactome obtained from NetworkAnalyst contains a total of 116 nodes with 341 interactions, where 20 nodes were from our common DEGs. The network comprises circular and hexagonal nodes representing TFs and common DEGs, respectively. Out of total 116 nodes, 111 nodes showed at least 3 TF-gene interactions. Notably, ICAM1, C3 and TYMP evidenced 49, 45 and 42 TF-gene interactions, respectively. The high number of interactions between the TFs and common DEGs are shown in Fig. 7 -A. Further network analysis with degree method revealed 30 most significant nodes including 15 TFs (IRF1, RERE, FOSL1, KLF8, MAZ, ZNF580, DRAP1, TRIM22, KLF16, WRNIP1, TFDP1, KLF11, CREB3L1, ZBTB11, and ZNF197) and their partner DEGs (Fig. 7-B). Similarly, the constructed gene-miRNA interaction network identified a total of 535 interactions among 121 nodes where 96 nodes represent miRNAs as shown in Fig. 8 -A. The gene ANLN exhibited maximum 59 interactions whereas BIRC3, PTGS2, INHBA, LIF, and ICAM1 were involved with 30 connections and/or more. All the nodes regulated at least three neighbours suggesting high level of miRNA association with the shared DEGs. Finally, topological study of the network identified 11 most important miRNAs (miR-1-3p, 34a-5p, 129-2-3p, 146a-5p, 16–5p, 101–3p, 671–5p, 155–5p, 124–3p, 20a-5p, and let-7b-5p) interacting with 19 DEGs (Fig. 8-B).

Fig. 7

TF-gene interaction network. The network was constructed using the shared DEGs and filtered with degree centrality. It shows (A) the interactions of 96 TFs with 20 DEGs, and (B) the 30 most significant nodes of the gene-TF network, which included 15 TFs and their interactions with 15 DEGs.The hexagonal (yellow) and circular (green) shaped nodes in the network indicate DEGs and TFs, respectively.

Fig. 8

Gene-miRNA interaction network. The network was constructed using the shared DEGs and filtered with degree centrality. It shows (A) the interactions of 96 miRNAs with 25 DEGs, and (B) the 30 most significant nodes of the gene-miRNA network that includes 11 miRNAs and their interactions with 19 DEGs.The hexagonal (yellow) and circular (blue) shaped nodes in the network indicate DEGs and miRNAs, respectively.

Candidate drug identification

Based on the shared DEGs between SARS-CoV-2 and T2D, we identified some potential drug molecules. We obtained 438 significant drugs/compounds to be associated with the DEGs that resulted from manual curation by considering FDR adjusted -value. Among them, estradiol and tetradioxin were the top 2 molecules based on the number of DEGs involved. Again, phencyclidine and profenamine molecules were found to be the most significant ones. The top 20 significant drug molecules along with the associated DEGs are provided in Table 1 .

Table 1

Top 20 significant drug candidates identified for shared DEGs between SARS-CoV-2 and T2D.

Drug/small molecule	p-value	Adj. p-value	Associated genes
Phencyclidine	1.16E-15	4.68E-12	IL1A, CCL20, IL1B, TNFAIP2, LIF, CXCL1, CXCL3, PTGS2, ICAM1, TYMP, BIRC3
Profenamine	1.94E-14	3.91E-11	IL1B, LIF, CXCL1, INHBA, CXCL3, PTGS2, ICAM1, BIRC3
8-Azaguanine	5.99E-14	8.04E-11	IL1A, CCL20, IL1B, LIF, CXCL1, INHBA, CXCL3, PTGS2, KRT6B, ICAM1, BIRC3
Nickel chloride	6.02E-13	6.06E-10	IL1A, CCL20, IL1B, TNFAIP2, CXCL1, INHBA, LTB, CXCL3, PTGS2, ICAM1, BIRC3
Chloropyramine	2.62E-12	2.11E-09	LIF, CXCL1, INHBA, PTGS2, ICAM1, BIRC3
Pizotifen	2.36E-11	1.58E-08	CXCL1, INHBA, CXCL3, PTGS2, ICAM1, BIRC3
MS-275	4.76E-11	2.74E-08	C3, IL1A, CCL20, IL1B, TNFAIP2, CXCL1, PTGS2, KRT6B, S100A8, BIRC3
0297417-0002B	1.53E-10	7.74E-08	LIF, CXCL1, LTB, PTGS2, KRT6B, ICAM1, BIRC3
Silica	1.89E-10	8.45E-08	ANLN, IL1A, EGR3, CCL20, IL1B, TNFAIP2, CXCL1, LTB, CXCL3, PTGS2, ICAM1, BIRC3
Gemcitabine	2.93E-10	1.07E-07	IL1A, CCL20, IL1B, CXCL1, CXCL3, ICAM1, TYMP, BIRC3
Peptidoglycan	3.313E-10	1.11E-07	IL1A, CCL20, IL1B, PTGS2, ICAM1
1-Nitropyrene	2.83E-10	1.14E-07	IL1A, CCL20, IL1B, CXCL1, LTB, CXCL3, ICAM1, BIRC3
Estradiol	8.55E-10	2.45E-07	EGR3, MME, CCL20, SLC6A14, ADAM22, LIF, KRT23, INHBA, CXCL3, PTGS2, ICAM1, C3, ANLN, IL1A, VNN1, IL1B, LTB, KRT6B, S100A8, BIRC3
Nickel sulphate	8.04E-10	2.49E-07	IL1A, CCL20, IL1B, CXCL1, INHBA, CXCL3, PTGS2, ICAM1, TYMP, BIRC3
3-Nitrofluoranthene	9.33E-10	2.50E-07	CCL20, CXCL1, LTB, CXCL3
Dexamethasone	1.12E-09	2.67E-07	IL1B, SLC6A14, LIF, LTB, PTGS2, S100A8, ICAM1, TYMP, BIRC3
Niclosamide	1.10E-09	2.77E-07	CCL20, CXCL1, INHBA, CXCL3, PTGS2, ICAM1, BIRC3
MG-132	1.35E-09	3.02E-07	IL1A, IL1B, CXCL1, PTGS2, ICAM1, TYMP, BIRC3
Simvastatin	2.24E-09	4.75E-07	IL1A, IL1B, LIF, LTB, PTGS2, S100A8, ICAM1, BIRC3
Thioridazine	3.22E-09	6.49E-07	LIF, INHBA, LTB, CXCL3, PTGS2, ICAM1, BIRC3

Discussion

This research work aimed to uncover the obscure genetic interactions between SARS-CoV-2 and T2D by maneuvering a series of bioinformatics approaches. To our best knowledge, this is the first endeavor to divulge the genomic and transcriptomic connections shared by the disease pair. In this study, we first performed the gene expression analysis of the SARS-CoV-2 and T2D in lung epithelium and pancreatic islet cells, respectively, as well as in peripheral blood using four different datasets. Considering the expression patterns, we identified 25 unique DEGs shared by the disease pair that manifested their genetic proximity. The subsequent course of actions in this study revolved around these shared DEGs employing functional enrichments in terms of GO, molecular pathways, protein-protein interactome, regulatory network, and drug target identification.

Functional enrichment analysis using the overlapping DEGs derived significant GO terms associated with diseases under consideration. Markedly, inflammation related signalling pathways, including cytokine activity, cytokine-mediated signalling pathway, chemokine activity, and cellular response to cytokine stimulus were found to be prevalent. Evidently, the comprehensive activity of these pathways may contribute to the cytokine storm, thus, playing a crucial role in the progression of COVID-19 [55]. A recent study revealed chronic inflammation as a prominent feature of T2D that may lead to abnormal blood clot formation, which explains the severity of T2D upon SARS-CoV-2 infection [56]. Another biological process is the regulation of cell proliferation, which is plausible since cell cycle regulation actively occurs in T2D patients to counteract insulin resistance by increasing insulin-secreting -cells despite the lack of success [57]. Therefore, the DEGs associated with this GO term may directly or indirectly involve in the regulation of -cell proliferation and hence, can be potential therapeutic targets [58].

In our GDA analysis, we found that most of the shared DEGs (18) were associated with breast carcinoma and cancer metastasis. The ICAM1 is crucial for infectious disease in viral replication modulation and has been evidenced to have key involvement in lung disease of severe COVID-19 patients [59]. Again, genetic and biological evidence suggest that ICAM1 has a vital role in the development of diabetes [60]. The dysregulation of ANLN (anillin actin binding protein) gene is associated with pulmonary carcinogenesis and it can be characterized as a prognostic biomarker and therapeutic candidate for lung cancer [61]. A recent study identified ANLN as a potential biomarker for bladder urothelial carcinoma [62]. Besides, 10 DEGs from this study were found to be linked to bladder cancers. A recent review highlighted the higher risk of COVID-19 and/or disease severity in individuals with breast cancer [63]. On the other hand, the least associated diseases were atherosclerosis, IBD, bladder and cervix cancers, and ulcerative colitis, which involved with 10 shared DEGs each. In particular, the intense inflammations due to IL-6, TNF-, and IL-1 can lead to aggravated atherosclerosis which in turn contribute to developing cardiovascular diseases [64]. Although no association was found between IBD and increased susceptibility towards SARS-CoV-2, a recent study suggested gut inflammation may contribute to viral entry by up-regulating ACE2 expression [65]. Other notable diseases were lung carcinoma, pancreatic carcinoma, rheumatoid arthritis, and asthma, which were also evident in some previous studies [[66], [67], [68], [69]]. Therefore, the associated DEGs could be of interest for risk management and genetic manipulation for therapeutic purposes.

Pathway analysis focusing on the common DEGs was carried out and identified similar pathways for COVID-19 and T2D. Amongst all, TNF signalling pathway, cytokine-cytokine receptor interaction, IL-17 signalling pathway, and photodynamic therapy-induced NF–B survival signalling were the top annotated pathways. Ouyang et al. found TNF-signalling pathway to be significantly enriched in severe COVID-19 conditions [70]. Recently, Akash et al. reviewed the direct involvement of TNF- with the pathogenesis of T2D. Activation of TNF- reduces the expression of glucose transporter type 4 (GLUT4) which is related to insulin metabolism. Besides, TNF- expression coincides with the inhibition of insulin receptor through serine phosphorylation of insulin receptor substrate-1 (IRS-1) [71]. Likewise, IL-17 pathway is involved in cytokine storm that regulates normal immune response [72]. Intriguingly, IL-17 also promotes insulin resistance by inducing angiotensin II type 1 receptor (AT1R). This lead to the production of nitric oxide (NO) in diabetic nephropathy resulting in complications in T2D patients [73]. Therefore, the increased disease severity of COVID-19 in T2D patients is not surprising. Now that we know the genetics of this complex interplay, blocking the involved DEGs could be an effective strategy in alleviating disease severity by regulating these signalling pathways. For instance, there are some clinical trials of antibodies against IL-17 for chronic inflammatory diseases such as psoriasis and arthritis [74]. These approaches can also be adopted to ameliorate the hyper-inflammation in SARS-CoV-2-infected T2D. Moreover, AT1R antagonists or vitamin D or in combinations are suggested to reduce the COVID-19 complications [75].

The PPI network analysis using the shared DEGs promisingly identified BIRC3, C3, MME, and IL1B as the potential hub-proteins, the proteins exhibiting a high degree of interactions and identified in multiple networks analyses. Notably, BIRC3 (Baculoviral IAP repeat containing 3) is connected with Marginal Zone B-Cell Lymphoma and is reported as NK cell immune checkpoint in tumor [76]. Likewise, C3 (Complement C3) is a principal element of the protective immune system against pathogens [77]. A recent study noticed hyper-activation of C3 in COVID-19 and suggested its suppression as a therapeutic strategy [78]. Importantly, C3 is suggested as a potential therapeutic target for diabetic nephropathy in a recent study [79]. Another hub-protein MME (Membrane metalloendopeptidase) in cerebrospinal fluid (CSF) was found to be involved in Alzheimer's disease progression [80]. Interestingly, in this study, we also identified Alzheimer's disease as the comorbidity of COVID-19. This protein directly employs neutrophil recruitment and consequently intensifies the immunopathology of COVID-19 [81]. The IL1B (Interleukin 1 beta) belongs to the interleukin 1 cytokine family and plays a vital role in arbitrating the inflammatory response. It is also involved in osteoarthritis pathogenesis [82], which was also evident in our comorbidity analysis. Moreover, 1L1B like pro-inflammatory cytokine level upraises in the bronchial alveolar lavage fluid of serious COVID-19 patients [83]. Consequently, the inflammatory response promoted by IL1B induces lung damage in COVID-19 [84]. Additionally, PTGS2 (Prostaglandin-endoperoxide synthase 2), exclusively detected by MCC method, is demonstrated to be a prognostic biomarker for gastric cancer [85]. Moreover, PTGS2 is lately reported as a key target molecule of the Qing-Fei-Pai-Du decoction (QFPDD), the traditional Chinese medicine that was successfully applied to treat COVID-19 [86]. Therefore, the identified hub-proteins could be considered as potential biomarkers for COVID-19. Further investigation could be conducted to confirm their biological involvement to highlight them as drug targets.

Furthermore, the regulatory network analyses elucidated notably high occurrence of TF-gene and gene-miRNA interactions. The ICAM1 (Intercellular adhesion molecule 1) gene was found to be involved with the highest number (49) of TF-gene interactions. In this study, we found interferon regulatory factor 1 (IRF1) as a potential TF. A previous study showed IRF1 directly regulates gene expression in response to herpes simplex virus (HSV) and MERS-CoV infections [87], probably by promoting IRF3 activation [88]. It is also responsible for the upregulation of interleukin-1 (IL1B) as we evidenced in our gene-TF network [89]. Also, krüppel-like factors (KLFs) are TFs that regulate various pathways and rendered metabolic abnormalities if altered. Notably, phenotypic expression related to KLF-11 includes decreased circulating insulin [90]. We identified KLF8, KLF11 and KLF16 as highly active TFs. Besides, we detected TF protein CREB3L1 which is in line with a previous study [91]. This TF is associated with ER stress [92] and PCSK1 gene, which was found to be linked with diabetes mellitus [93]. Another TF protein FOSL1 interacts with CCL20 gene, which was found upregulated in this study. The elevated level of CCL20 is implicated in cytokine storm of COVID-19, obesity, and insulin resistance [94,95]. Also, the CCL20 is regulated by the NF–B subunits in pancreatic -cells [95]. Importantly, FOSL1 is reported to negatively regulate the type I interferon (IFN–I) response in host against malaria and viral infections and is identified as a potential drug target for controlling malaria and other diseases [96].

We also detected 11 potential miRNAs that might have shared pathogenetic relationship between COVID-19 and diabetes. For example, miR-1-3p and 20a-5p identified in this study reported to be involved in viral respiratory diseases and can be used to design anti-viral drugs [97]. Another miRNA miR-34a-5p increases high glucose-mediated apoptosis in cardiomyocytes by reducing anti-apoptotic BCL2 protein [98]. Thus, the upregulation of miR-34a-5p in the diabetic condition is most likely to increase apoptosis [99]. Interestingly, the elevated level of miRNA-34a was also observed in patients who experienced acute myocardial infarction (MI) leading to heart failure [99]. Further, some miRNAs associated with lung diseases (e.g., let-7b-5p) and asthma (e.g., miR-155–5p, miR-16–5p) were highly prevalent in COVID-19 [100]. Furthermore, miR-16 and miR-155 modulate inflammatory responses by targeting the MCL1 and NFKBIA hub-proteins [100,101]. Therefore, the association of identified TFs and miRNAs with COVID-19 pathogenesis and diabetes could reveal other potential drug targets with further study.

Moreover, taking the common DEGs into consideration with the DSigDB database, several drug/molecule candidates were proposed. As expected, the top interacted component estradiol was found to have a positive protective effect on COVID-19 patients in several experimental investigations that suggests its prospect as a therapeutic candidate against SARS-CoV-2 [102]. Again, dexamethasone which had association with 9 shared DEGs, was found effective in COVID-19 patients who require oxygen supplement and/or mechanical ventilation [103]. Another drug candidate tetradioxin (not shown in Table 1) showed interactions with 18 shared DEGs. Previously, dioxin showed immune suppression activity through the inhibition of CD4+ T-cell differentiation into Th1, Th2 and Th17 effector cells [104]. Strikingly, activation of Th17 is mainly responsible for IL-17 release, the prominent cause of cytokine storm in COVID-19 [105,106]. Therefore, the identified drug candidates can further be assessed for their potency in alleviating the COVID-19 complication. Altogether, the obtained results demand clinical expeditions that may improve our understanding of the genetic impact of T2D on SARS-CoV-2 as well as to propose putative therapeutic targets.

Conclusions

To conclude, the pathophysiological outcome of SARS-CoV-2 infection significantly differs in different people based on their underlying health conditions. In this study, we explored the molecular pathogenesis of COVID-19 in patients with pre-existing T2D hoping to find the genetic determinants that cause higher infection risk and disease severity. In doing so, we observed that both diseases exhibited similarity in gene expression profile to some extent. We also identified several key proteins that could be the source of potential biomarkers in early disease prognosis. In addition, the identified regulatory molecules such as TFs and miRNAs may highlight the disease-modifying factors involved in the COVID-19 progression. These regulatory molecules along with the signalling pathways can be targeted for the development of novel COVID-19 therapeutics. Furthermore, the identified potential drug candidates need experimental validation for their therapeutic evaluation. Thus, we believe that this study will shed light on the understanding of the intertwined pathobiology of the complexity of SARS-CoV-2 infected T2D patients.

Contributions

All authors contributed to the manuscript. M.B.I., U.N.C., and M.A.M. conceived and designed the study. M.B.I. and U.N.C. analyzed the data and wrote the R programming code for the development of the pipeline. M.B.I., U.N.C., Z.N. and M.A.M. conducted all other bioinformatic analyses and wrote the manuscript. S.U. and M.B.A. provided insightful guidance and feedback for the study. All authors read and approved the final version of the manuscript.

Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.