Introduction
Type 2 diabetes (T2D) has rapidly emerged as a global health issue, driven by factors such as population growth, aging, urbanization, and the rising prevalence of obesity and physical inactivity [1]. The International Diabetes Federation (IDF) estimated that globally in 2021 10.5% of individuals (536.6 million people) aged 20–79 years had diabetes, a figure that is expected to increase to 12.2% (783.2 million) by 2045. Globally more than half a billion individuals live with diabetes, accounting for nearly 10.5% of the world adult population [2].
Disease comorbidity refers to the coexistence of one or more additional diseases or conditions alongside a primary or index disease. [3]. Long-term presence of T2D increases the risk of developing comorbidities such as obesity, hypertension, diabetic nephropathy (DN), endothelial dysfunction, end-stage renal disease, etc., which can impede antihyperglycemic therapies and overall management of T2D [4, 5]. A significant proportion of individuals with diabetes are obese (approximately 90%), while between 20% and 50% of individuals who have diabetes develop DN [6, 7]. Current IDF statistics suggest that every 6 seconds, a person dies from diabetes or its complications, with 50% of these deaths (totalling 4 million annually) occurring in individuals under 60 years old [8].
T2D is a metabolic disorder defined by insulin resistance and dysfunction of pancreatic β-cells, resulting in chronic hyperglycemia [9]. Chronic hyperglycemia in T2D is linked to damage, dysfunction, and failure of various organs, including the retina, kidneys, nervous system, heart, and blood vessels [10]. Both T2D and obesity are linked to insulin resistance [10]. DN is among the most common and serious complications of T2D and is linked to higher rates of morbidity and mortality in individuals with diabetes [11]. Clinical studies, including epidemiological analysis, suggest pathogenic links between T2D and obesity as well as between T2D and DN. In a study by Ruze et al. (2023) [12], it has been reported that obesity, characterized by shared genetic and environmental factors in its development, exacerbates the influence of genetic susceptibility and environmental factors on the pathophysiology of T2D. Raman et al. (2010) [13] reported the prevalence of metabolic syndrome, characterized by central obesity, glucose intolerance, high insulin levels, low HDL cholesterol, elevated triglycerides, and hypertension, and its impact on microvascular complications such as diabetic nephropathy, diabetic retinopathy, and diabetic neuropathy in the Indian population with T2D [13].
Complex diseases arise from a combination of factors, often accompanied by genetic perturbations. These diseases can be explored using a network systems biology approach because network biology has established itself as a potent methodology for the analysis of complex diseases, focusing on studying networks of co-expressed genes and their associated functional modules. Network-centric approaches for studying complex diseases hold significant promise in identifying critical targets that function as both biomarkers and therapeutic targets [14, 15]. Numerous transcriptomic and genetic studies have been conducted on T2D. CD44 and CCL5 genes have been shown to be associated with wound healing in diabetic patients through analysis of transcriptome of dermal lymphatic endothelial cells [16]. A genetic study by Ali (2013) [17] reported several genes, such as PPARG, IRS1, IRS-2, KCNJ11, WFS1, HNF1A, HNF1B, and HNF4A, to be associated with T2D that were involved in insulin secretion and action. A transcriptomic study by Tonyan et al. (2022) [18] identified various differentially expressed genes (DEGs) involved in pathways such as lipid metabolism, immune response, and the ubiquitin-proteasome system in T2D.
Several studies have demonstrated a direct relationship between T2D and obesity, as well as T2D and DN. Bima et al. (2022) [19] identified the ERBB2, FN1, FYN, HSPA1A, HBA1, and ITGB1 genes as potential genetic markers linking T2D and obesity using a computational network biology approach. In another study, the PTPRC, CD53, IRF8, IL10RA, and LAPTM5 genes were found to play essential roles in the molecular mechanisms of DN, a complication of T2D [20].
An epidemiological study to identify diabetic nephropathy risk factors in T2D obese people based on community T2D patients reported an association between the 3 clinical conditions [21]. However, most genetic, transcriptomics, and bioinformatics studies have focused on the direct association of T2D with obesity or T2D with DN, and there are no reports on the relationship between T2D, DN, and obesity at the molecular level. This leaves a gap in the understanding of the molecular underpinning genes and pathways common in co-occurrence. The present study was therefore undertaken to identify the association at the genetic level between T2D, obesity, and DN using an in silico network systems biology approach.
Materials and methods
Data retrieval
To investigate the association of T2D with DN and obesity, the gene expression profile of 3 different microarray datasets were retrieved from the NCBI Gene Expression Omnibus (GEO) database using the GEO query package of R programming language [22]. The T2D dataset with accession number GSE20966 [23] was on account of GPL1352 [U133_X3P] Affymetrix Human X3P Array platform including beta cells from the pancreas of 10 T2D patients and 10 healthy patients. The DN microarray dataset GSE1009 [24], comprising 3 diabetic nephropathy patients and 3 control subjects, was based on the GPL8300 [HG_U95Av2] Affymetrix Human Genome U95 Version 2 Array platform and provided gene expression profiles of human kidney glomeruli. The obesity dataset, GSE9624 [25], was a GPL570 [HG-U133_Plus_2] Affymetrix Human Genome U133 Plus 2.0 Array gene expression profile of adipose tissue from 5 obese and 6 normal weight patients.
Pre-processing of microarray data and DEG identification
The 3 GEO datasets underwent pre-processing steps that included data consolidation, normalization via log2 transformation, and conversion of probe IDs to gene symbols using the gprofiler2 package in RStudio [26]. The data were corrected and normalized, followed by the generation of a box plot. Volcano plots for all 3 datasets were constructed and visualized using the ggplot2 package [27]. To identify differential gene expression between disease and control states, up-regulated genes were obtained with p-value < 0.05 and logFC > 1, while down-regulated genes were obtained with p-value < 0.05 and logFC < –1, utilizing the limma package of R [28].
The 3 datasets of DEGs obtained for T2D, obesity, and DN were analysed using InteractiVenn software, an online tool that generates Venn diagrams to visualize overlaps and intersections among multiple datasets [29], which was used to identify commonly expressed DEGs amongst the 3 comorbidites.
PPI network construction, analysis, and gene enrichment
STRING (Search Tool for the Retrieval of Interacting Genes/Proteins) web-based visualization software, was used to construct and analyse the protein-protein interaction (PPI) network [30], which was further analyzed using Cytoscape [31]. Topological analysis and identification of functional modules of the networks was performed using the Cytoscape plug-ins CytoHubba [32] and CytoCluster [33].
Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis was performed using Enrichr-KG (Enrichr Knowledge Graph), which returns an integrated network of the top enriched terms from multiple libraries, connected to their overlapping genes [34].
Results
Data retrieval and identification of DEGs
The microarray dataset for T2D, obesity, and DN with series identifiers GSE20966, GSE9624, and GSE1009, respectively, were retrieved from the NCBI GEO database and pre-processed to identify the DEGs. Data pre-processing steps involved data consolidation and normalization along with selection and filtering the sample data using R studio. The volcano plot for all the three datasets were generated using ggplot2 package of R studio (Fig. 1).
To identify the differential gene expression between disease and control states, p-value < 0.05 and logFC > 1 was used to obtain the up-regulated genes while p-value < 0.05 and logFC < –1 was used to obtain the down-regulated genes using the limma package. In T2D, a total of 668 DEGs were identified in which 443 were up-regulated and 225 were down-regulated. Compared to normal patients, 2671 genes were significant DEGs in obese patients with 1249 up-regulated genes and 1422 down-regulated genes. In the DN patient group, a total of 2471 genes were significant DEGs, of which 1643 were up-regulated and 828 were down-regulated. A cross-comparison analysis of InteractiVenn identified the common DEGs amongst T2D, DN, and obesity, as shown in Figure 2, which resulted in the identification of 13 common DEGs.
PPI construction, analysis, and identification of hub genes
The 13 common DEGs shared amongst T2D, DN and obesity were imported in STRING database to construct PPI network. The network obtained from STRING contained 93 nodes (representing genes) and 866 edges (representing interaction between nodes), as shown in Figure 3.
The top 10 hub genes were identified and ranked according to their significance, i.e., highly significant genes in red, followed by orange, and least significant genes in yellow using the CytoHubba (Maximum Clique Centrality method) Cytoscape plug-in (Fig. 4A). The highly significant hub genes identified were AKR1C3, CYP19A1, AKR1D1, HSD17B3, and HSD17B6 followed by their validation using the ClusterONE algorithm of the Cytoscape plug-in CytoCluster. Four clusters were obtained (Tab. 1) out of which 3 clusters with p-value < 0.05 were merged (Fig. 5) and further analysed using CytoHubba to obtain the top 10 hub genes, namely AKR1C3, CYP19A1, AKR1D1, HSD17B3, HSD17B8, and HSD17B2 (Fig. 4B).
After a cross comparison between the highly significant genes obtained using Cytoscape plug-ins CytoHubba and CytoCluster, 4 hub genes were identified, i.e., AKR1C3, CYP19A1, AKR1D1, and HSD17B3.
Gene set enrichment analysis
The gene products AKR1D1, HSD17B3, AKR1C3, and CYP19A1 were found to be members of the KEGG pathway of steroid hormone biosynthesis. AKR1C3 was also part of ovarian steroidogenesis, arachidonic acid metabolism, and folate biosynthesis, while CYP19A1 and AKR1D1 were members of ovarian steroidogenesis and primary bile acid biosynthesis, respectively (Fig. 6).
Discussion
In the present study, 4 hub genes were identified as underlying molecular factors responsible for T2D, obesity, and DN by analyzing DEGs associated with these 3 clinical conditions. The hub genes identified in the study were further investigated using literature mining to validate their association with the pathways in which the hub genes were enriched.
AKR1C3 (aldo-keto reductase 1C3) identified in the present study has been reported to be responsible for producing potent androgens in peripheral tissues [35] and plays a crucial role in steroid metabolism by converting the precursor androstenedione into testosterone [36, 37]. Franko et al. (2020) reported that hyperglycemia increases the mRNA levels of AKR1C3 in T2D [38] while Svensson et al. (2008) reported the involvement of AKR1C3 with obesity [39].
CYP19A1 (cytochrome P450, family 19, subfamily A, peptide 1) has been reported by Baddela et al. (2020) to be involved in activating p-Akt/HIF-1α signaling pathway, which in turn regulates the metabolism of glucose and lipids and is therefore associated with T2D [40]. AKR1D1 (aldo-keto reductase 1D1) has been reported by Nikolaou et al. (2019) [41] to play a crucial role in bile acid synthesis, and its silencing enhances insulin sensitivity and leads to lipid accumulation, thereby promoting hepatocyte inflammation. HSD17Bs (17β-hydroxysteroid dehydrogenases) comprises a group of 15 enzymes, and most members of this family, including HSD17B3, are involved in regulating the biological activity of sex hormones [42]. It was reported by Wang et al. (2023) [43] that HSD17B2 and HSD17B3 are located in the endoplasmic reticulum and are involved in lipid biosynthesis, lipid metabolism, and steroid biosynthesis.
Sex steroid hormones significantly influence metabolic pathways, playing a pivotal role in the development of metabolic disorders, including T2D. According to Zhang et al. (2024) [44], these hormones affect metabolism, adipose tissue distribution, and expansion by binding to specific receptors within the adipose tissue along with a reduction in estrogen and/or androgens commonly associated with central obesity.
The present study explored the molecular association between T2D, obesity, and DN using a network systems biology approach. This study identified 4 hub genes, namely AKR1C3, CYP19A1, AKR1D1, and HSD17B3, predominantly involved in steroid hormone biosynthesis pathway, which may contribute to the disruption of metabolic pathways, ultimately leading to the comorbidities associated with T2D.
Our findings highlighted the common underlying genes and pathways amongst these diseases which may account for their comorbidity. Although our study aimed to map key genes and pathways between T2D, obesity, and DN, our results at present are observational, based on experimental datasets of microarray analysis. Therefore, a limitation of our study is that the identified genes have not been validated by clinical or experimental studies, which could provide valuable insights for the development of common therapeutic strategies for T2D, obesity, and DN.
Conclusions
This study is the first of its kind to investigate the molecular association between T2D, obesity, and DN in terms of genes and pathways. Genes common amongst the comorbidities, namely AKR1C3, CYP19A1, AKR1D1, and HSD17B3, were identified, which were found to be involved in the same pathway of steroid hormone biosynthesis. Further experimental validation and clinical studies need to be carried out and may become a focus of future research.
Article information
Availability of data and materials
The datasets generated or analyzed during this study are available from the corresponding author upon reasonable request.
Ethics approval and consent to participate
The study did not involve human subjects or animal models, and therefore no ethical clearance was required.
Author contributions
TRS: conceived and designed the study; M: collected the data, performed the analysis, and wrote the first draft; TRS interpreted the data and performed major revisions of the first draft. Both authors reviewed and edited the manuscript and approved the submission.
Funding
The authors did not receive any funding for this study.
Acknowledgements
The authors would like to thank their parent institute Panjab University Chandigarh, India for providing the infrastructure for carrying out the research.
Conflict of interest
The authors declare no conflict of interest.