Identification of Crucial Genes and Pathways Associated with Atherosclerotic Plaque in Diabetic Patients

Background Patients with diabetes have more calcification in atherosclerotic plaque and a higher occurrence of secondary cardiovascular events than patients without diabetes. The objective of this study was to identify crucial genes involved in the development of diabetic atherosclerotic plaque using a bioinformatics approach. Methods Microarray dataset GSE118481 was downloaded from the Gene Expression Omnibus (GEO) database; the dataset included 6 patients with diabetic atherosclerotic plaque (DBT) and 6 nondiabetic patients with atherosclerotic plaque (Ctrl). Differentially expressed genes (DEG) between the DBT and Ctrl groups were identified and then subjected to functional enrichment analysis. Based on the enriched pathways of DEGs, diabetic atherosclerotic plaque-related pathways were screened using the comparative toxicogenomics database (CTD). We then constructed a protein–protein interaction (PPI) network and transcription factor (TF)–miRNA–mRNA network. Results A total of 243 DEGs were obtained in the DBT group compared with the Ctrl group, including 85 up-regulated and 158 down-regulated DEGs. Functional enrichment analysis showed that up-regulated DEGs were mainly enriched in isoprenoid metabolic process, DNA-binding TF activity, and response to virus. Additionally, DEGs participating in the toll-like receptor signaling pathway were closely related to diabetes, carotid stenosis, and insulin resistance. The TF–miRNA–mRNA network showed that toll-like receptor 4 (TLR4), BCL2-like 11 (BCL2L11), and glutamate-cysteine ligase catalytic subunit (GCLC) were hub genes. Furthermore, TLR4 was regulated by TF signal transducer and activator of transcription 6 (STAT6); BCL2L11 was targeted by hsa-miR-24-3p; and GCLC was regulated by nuclear factor, erythroid 2 like 2 (NFE2L2). Conclusion Identification of hub genes and pathways increased our understanding of the molecular mechanisms underlying the atherosclerotic plaque in patients with or without diabetes. These crucial genes (TLR4, BC2L11, and GCLC) might function as molecular biomarkers for diabetic atherosclerotic plaque.


Introduction
Atherosclerosis is a chronic inflammation disease and the leading cause of morbidity and mortality globally. 1 Atherosclerosis is a slowly progressive process, characterized by an accumulation of lipid in the arterial wall accompanied by multifocal structural alterations, leading to atheromatous plaque formation. 2,3 Over time, the large necrotic lipid core is covered by a fibrous cap until, in advanced stages, the stability of the cap is destroyed, inducing plaque rupture and thrombosis, which can manifest as stroke or myocardial infarction. Cardiovascular disease is the leading cause of death in patients with diabetes. 4 There is increasing evidence that diabetes induces hypercoagulability, which has a role in plaque rupture and increases the incidence and severity of clinical events. 4 Previous researchers have reported the relationship between plaque characteristics and patients with and without diabetes. Burke et al indicated that total plaque in diabetic patients was significantly greater than that of nondiabetic individuals; in addition, the inflammatory response of diabetic plaques was stronger than that of nondiabetic plaque. 5 An optical coherence tomography imaging study by Kato et al revealed that plaques in diabetic patients had a higher incidence of calcification and thrombus. 6 Furthermore, van Haelst et al found that patients with diabetes had more calcification in atherosclerotic plaque and a higher occurrence of secondary cardiovascular events than patients without diabetes. 7 Even though we can distinguish atherosclerotic plaque in diabetic and nondiabetic patients from morphologic fields, the effect of diabetes on gene expression in atherosclerotic plaque is not fully understood.
Macrophage accumulation plays a vital role in both plaque progression and stability, which can promote inflammation and aggravate disease. 8,9 Thus, we selected a gene expression dataset (GSE118481) containing diabetic plaque macrophage (DBT) and nondiabetic plaque macrophage (Ctrl) for analysis. Differentially expressed genes (DEGs) between the DBT and Ctrl groups were identified, functional enrichment analysis of the DEGs was performed, and disease-related pathways were screened. We then constructed a protein-protein interaction (PPI) network and sub-network. Subsequently, microRNA (miRNA) and transcription factors (TFs) of DEGs were predicted and an integrated TF-miRNA-mRNA network was constructed. The analysis process of this study is shown in Supplementary Figure 1. We aimed to further understand the molecular mechanism by which diabetes promotes the formation of atherosclerotic plaque and to determine potential gene targets for personalized diagnosis and treatment strategies of diabetic atherosclerotic plaque.

Data Sources
The gene expression profile GSE118481 based on the GPL10558 Illumina HumanHT-12 V4.0 expression BeadChip platform was downloaded from the Gene Expression Omnibus (GEO) database (website: http:// www.ncbi.nlm.nih.gov/geo/). 10 This microarray data set included 16 nondiabetic samples (6 asymptomatic and 10 symptomatic) and 8 diabetic plaque samples (6 asymptomatic and 2 symptomatic). In order to study the effect of diabetes on atherosclerotic plaque, we selected asymptomatic samples for subsequent analysis. Therefore, 12 samples (6 DBT and 6 Ctrl) were included, and the clinical characteristics of these patients are listed in Supplementary Table 1. There were no significant differences in age (P = 0.24) and sex (P > 0.05) between the two groups.

Data Preprocessing and Identification of DEGs
The series matrix file for GSE118481 dataset was obtained from the GEO database, 10 and the expression data of 6 DBT and 6 Ctrl macrophage samples were extracted for further analysis. Microarray expression profiling was standardized by Bioconductor bead array package, 11 and the distribution of expression in each sample was visualized by boxplots. The probe ID was converted to a gene symbol using the annotation file, and probes that did not mapped to gene symbols were removed. If multiple probes matched the same gene, the mean value of probes was calculated. Empirical Bayes moderated t-test in the limma package (version 3.40.6) 12 was used to identify DEGs between DBT and Ctrl samples. DEGs with P < 0.05 and |log fold change (FC)| >0.585 were considered statistically significant. The ggplot2 and heatmap of R (http://www. R-project.org/) were utilized to visualize the DEGs.

Disease-Related Co-Expression Network Analysis
Using Pearson correlation coefficients (r) in the stats of R package (version 3.6.1; http://www.R-project.org/), the co-expression of DEGs in DBT and Ctrl samples was, respectively, analyzed. Pairs with r > 0.95 and P < 0.05 were selected for further study. Cytoscape was applied to construct a co-expression network of Ctrl and DBT groups, and then the topological properties of the two networks were analyzed by using CytoNCA in Cytoscape. 13 Furthermore, t-tests were used to calculate the difference between Ctrl and DBT networks. The sub-networks of DBT group were structured using the Cytoscape MCODE plugin, 14 and networks with a score >3 were selected.

Functional Enrichment Analysis of DEGs
To understand the major biological functions of DEGs, we analyzed Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways of up-and downregulated DEGs using the clusterProfiler package. Significant enrichment was defined by P < 0.05 and count >2.

Disease-Related Pathway Screening
Diabetes, carotid stenosis, and insulin resistance have effects on the development of diabetic atherosclerotic plaque. Therefore, to further identify diabetic plaque-related pathways, we screened pathways relevant to these diseases from the comparative toxicogenomics database (CTD), and these pathways were integrated with the KEGG pathways in the previous step.

PPI Network Construction
To determine the relationships among DEGs, we mapped the DEGs to the Search Tool for the Retrieval of Interacting Genes (STRING, version 11.0, http://www.string-db.org/) database, and only interactions with a combined score >0.7 were selected. Then, Cytoscape software was used to establish the PPI network, and hub nodes in the network were identified by CytoNCA. The PPI network was further analyzed by MCODE to explore functional modules, and score >5 was selected as the threshold.

TF-miRNA-mRNA Regulatory Network Analysis
To further understand the regulatory mechanism of DEGs, miRNAs of target genes were predicted using four available databases: miRWalk3.0, 15 TargetScan, 16 MiRDB, 17 and MirTarBase. 18 Putative miRNAs with score >0.95 and supported by at least two databases were selected; additionally, TF-target interactions were predicted by Transcriptional Regulatory Relationships Unraveled by Sentence-based Text mining (TRRUST) (https://www.grnpedia.org/trrust/). 19 Subsequently, miRNA-target pairs and TF-target pairs were integrated to construct the TF-miRNA-mRNA regulatory network.

Identification and Analysis of DEGs
The raw data were processed and the boxplots showed good normalized properties ( Figure 1A). A total of 33,984 probes were obtained after annotation. A total of 243 DEGs were identified between the DBT and Ctrl groups, including 85 up-regulated and 158 downregulated DEGs. The heat map and volcano plot of DEGs are shown in Figure 1B and C.

Disease-Related Co-Expression Network Analysis
Co-expression network analysis showed that 144 pairs and 162 DEGs were identified in the Ctrl group (Figure 2A), and 191 relationships and 170 DEGs were screened in the DBT group ( Figure 2B). The topological properties of the Ctrl and DBT networks indicated that betweenness, closeness, and degree values were significantly higher in the DBT group than in the Ctrl group ( Figure 2C), suggesting that co-expression of DEGs was more abundant in the DBT group than in the Ctrl group. The sub-network (score = 3.333) of the DBT group was constructed and included serine/threonine kinase 32B (STK32B), tripartite motif containing 22 (TRIM22), ninjurin 2 (NINJ2), and transmembrane protein 114 (TME114) ( Figure 2D).

Functional Enrichment Analysis
Functional enrichment analysis of up-regulated and downregulated DEGs was performed using the clusterProfiler tool. The top 10 significantly enriched GO terms and KEGG pathways are shown in Figure 3A and B. The up-regulated DEGs were significantly enriched in GO terms related to isoprenoid metabolic process, response to estradiol, and retinal metabolic process, and the down-regulated DEGs were markedly associated with regulation of DNA-binding TF activity, alphaamino metabolic process, and response to virus ( Figure 3A). For the KEGG pathway analysis, up-regulated DEGs were primarily involved in adipocytokine signaling pathway, cosphingolipid biosynthesis-lacto and neolacto series, and non-alcoholic fatty liver disease; in addition, downregulated DEGs were mainly involved in folate biosynthesis, cysteine and methionine metabolism, and Chagas disease (American trypanosomiasis) pathways ( Figure 3B).

Disease-Related Pathway Analysis
A total of 13 pathways were closely related to diabetes mellitus, one pathway was associated with carotid stenosis, and 10 pathways were associated with insulin resistance (Table 1). All three of these diseases were relevant to the toll-like receptor signaling pathway. Genes such as toll-like receptor 8 (TLR8), toll-like receptor 4 (TLR4), mitogen-activated protein kinase 4 (MAP2K4), and interferon regulatory factor 5 (IRF5) were involved in this pathway. In addition, 10 pathways were

TF-miRNA-mRNA Regulatory Network Analysis
After screening, 102 miRNA-mRNA pairs and 114 TF-mRNA pairs were predicted, and then these pairs were integrated to structure a TF-miRNA-mRNA regulatory network. A total of 154 interactions were identified, involving 57 genes, 30 miRNA, and 75 TFs ( Figure 5). In this regulatory network, we noted glutamate-cysteine ligase catalytic subunit (GCLC), BCL2 like 11 (BCL2L11), and TLR4 had higher degrees. GCLC was targeted by the TF nuclear factor, erythroid 2 like 2 (NFE2L2); BCL2L11 was targeted by hsa-miR-24-3p, and regulated by TF forkhead box O3 (FOXO3), and TLR4 was related to the process of three diseases and regulated by TF signal transducer and activator of transcription 6 (STAT6).

Discussion
Diabetes is known to be associated with atherosclerotic plaque; however, the underlying molecular mechanism of the effect of diabetes on atherosclerotic plaque has not been fully elucidated. We analyzed gene expression patterns involved in diabetic atherosclerotic plaque using the dataset GSE118481. The results revealed that the tolllike receptor signaling pathway was associated with the pathogenesis of diabetic atherosclerotic plaque. Additionally, TLR4, BCL2L11, and GCLC were potential biomarkers for atherosclerotic plaque in patients with diabetes.
The formation and progression of atherosclerotic plaque is related to the accumulation of monocyte-derived macrophage in the arterial wall. Compared with patients without diabetes, the plaques in the coronary arteries of patients with diabetes gene generally exhibit larger necrotic cores and significantly greater inflammation, mainly composed of macrophages and T lymphocytes. 20 Based on analysis of disease-related pathways, we observed that the toll-like receptor signaling pathway was significantly associated with diabetes mellitus, carotid stenosis, and insulin resistance; additionally, hub gene TLR4 was involved in this pathway. Madhur et al 21 indicated that inflammation response could reduce the stability of atherosclerotic plaques in animal models. It is reported that the TLR signaling pathway is associated with systemic inflammation and immune response, and participates in angiogenesis, survival, and repair. 22,23 Meanwhile, TLR4, as a member of the TLR family, is believed to activate nuclear factor-κB in response to short-chain fatty acids, triggering further activation of the immune system. 24 Thus, TLR4 induced inflammation plays an important role in atherosclerotic plaque stability. Xu et al demonstrated that TLR4 was preferentially expressed by macrophages in human lipid-rich atherosclerotic lesions, where it might play a role to enhance and maintain the innate immunity and inflammation. In addition, the up-regulation of TLR4 in macrophages by oxidized low-density lipoprotein (LDL) suggested that TLR4 might provide a potential pathophysiological link between lipids as well as inflammation and atherosclerosis. 25 In this analysis, we also found the relationship between TLR4 and diabetes. Devaraj et al demonstrated that expression of TLR4 was significantly increased in patients with type 1 diabetes, suggesting that TLR4 contributes to the proinflammatory state in diabetes. 26 Moreover, knockout of TLR4 might alleviate inflammation in diabetic rats 27 and TLR4 antagonist could attenuate atherogenesis in mice with diabetes. 28 The antidiabetic drug class thiazolidinediones (TZDs) have been reported to reduce the risk of atherosclerosis in patients with type 2 diabetes, which might have an anti-atherosclerotic effect by inhibiting the TLR4 signaling pathway. 29 These findings emphasized the importance of TLR4 in plaque formation of patients with diabetes. In the present study, we found that TLR4 was regulated by the TF

DovePress
STAT6. STAT6 has a major role in the immune system 30 and is associated with macrophage polarization, which is critically involved in atherosclerosis progression and regression. 31 Based on our results, we speculated that TLR4 and STAT6 might participate in the pathogenesis of diabetic atherosclerotic plaque via the TLR signaling pathway.
In the regulatory network, BCL2L11 had higher degree and was considered a hub gene. BCL2L11 encodes BCL-2 protein family, and its members participate in various cellular activities as anti-or proapoptotic regulators. 32 A previous study revealed that BCL2L11 was connected to apoptosis of podocytes in diabetes. 33 However, there are few reports about the

216
association between BCL2L11 and carotid plaque. Our analysis showed that BCL2L11 was targeted by hsa-miR-24-3p. Erener et al observed that miR-24-3p was an effective biomarker to predict and diagnose diabetes. 34 Moreover, miR-24-3p was found to limit macrophage vascular inflammation and slow the progression of atherosclerotic plaque. 35 Taken together, hsa-miR-24-3p might affect the progression of diabetic atherosclerotic plaque by directly targeting BCL2L11. However, the specific regulatory mechanism of BCL2L11 in diabetic atherosclerotic plaque needs further elaboration.
We also found that GCLC was closely related to diabetic atherosclerotic plaque. GCLC is a rate-limiting enzyme of glutathione synthesis, and it is involved in susceptibility to myocardial infarction. 36 Callegari et al found that the gain and loss of the ability to synthesize glutathione especially in macrophages had reciprocal effects on the initiation and progression of atherosclerosis at multiple sites in apoE-/-mice. 37 Jain et al reported that the plasma level of GCLC was lower in diabetic patients than in healthy controls. 38 Moreover, the GCLC polymorphism was associated with cellular redox imbalances and modulate the risk for diabetic nephropathy. 39 In the

217
TF-miRNA-mRNA network, GCLC was regulated by NFE2L2 (also known as NRF2), which is considered to be a master regulator of the antioxidant response. 40 It regulates the expression of several genes including Phase II metabolic and antioxidant enzymes, and therefore plays an important role in preventing oxidative stress-mediated diabetes and related complications. 41 In addition, overexpression of Nrf2 could protect pancreatic cells from oxidative damage in diabetes. 42 Furthermore, a study by Harada et al revealed that activation of Nrf2 was observed in advanced atherosclerotic plaques, suggesting that Nrf2 might influence the inflammatory reactions in the plaques. 43 Thus, we speculated that GCLC targeted by NFE2L2 might participate in the pathogenesis of diabetic atherosclerotic plaque.
Some limitations should be noted in the current study. First, the sample size of this study was small; further investigations based on a larger sample should be performed. Second, hub genes were identified using bioinformatics analysis; thus, experimental studies are needed to validate our results. Despite these limitations, this study provided some new insights into the pathogenesis and treatment of diabetic atherosclerotic plaques. Further large-scale studies are needed to corroborate these findings and investigate the potential underlying mechanisms involved. Meanwhile, clinical trials with more detailed investigation are also warranted before genes such as TLR4, BCL2L11, GCLC can be used in clinical setting.

Conclusion
In summary, we have conducted a comprehensive bioinformatics analysis of DEGs between diabetic and nondiabetic atherosclerotic plaque. Several genes have been identified with different expression patterns in diabetic and non-diabetic atherosclerotic plaque, such as TLR4, BCL2L11, GCLC, STAT6, and NFE2L2, as well as hsa-miR-24-3p. Meanwhile, pathway analysis showed that these genes were involved in the toll-like receptor signaling pathway. These findings provided better understanding of the underlying molecular mechanisms of diabetic atherosclerotic plaque. However, further research of these candidate genes was needed to confirm their effects in diabetic atherosclerotic plaque.

Publish your work in this journal
Pharmacogenomics and Personalized Medicine is an international, peer-reviewed, open access journal characterizing the influence of genotype on pharmacology leading to the development of personalized treatment programs and individualized drug selection for improved safety, efficacy and sustainability. This journal is indexed on the American Chemical Society's Chemical Abstracts Service (CAS). The manuscript management system is completely online and includes a very quick and fair peer-review system, which is all easy to use. Visit http://www.dovepress.com/testimonials.php to read real quotes from published authors.