Integrated Bioinformatics-Based Identification of Potential Diagnostic Biomarkers Associated with Diabetic Foot Ulcer Development.

The present study was designed to detect possible biomarkers associated with diabetic foot ulcer (DFU) incidence in an effort to develop novel treatments for this condition. The GSE7014 and GSE29221 gene expression datasets were downloaded from the Gene Expression Omnibus (GEO) database, after which differentially expressed genes (DEGs) were identified between DFU and healthy samples. These DEGs were then arranged into a protein-protein interaction (PPI) network, and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway and Gene Ontology (GO) term enrichment analyses were performed to explore the functional roles of these genes. In total, 1192 DEGs were identified in the GSE7014 dataset (900 upregulated, 292 downregulated), while 1177 were identified in the GSE29221 dataset (257 upregulated, 919 downregulated). GO analyses revealed these DEGs to be significantly enriched in biological processes including sarcomere organization, muscle filament sliding, and the regulation of cardiac conduction, molecular functions including structural constituent of muscle, protein binding, and calcium ion binding, and cellular components including Z disc, myosin filament, and M band. These DEGs were also enriched in the adrenergic signaling in cardiomyoctes, dilated cardiomyopathy, and tight junction KEGG pathways. Together, the findings of these bioinformatics analyses thus identified key hub genes associated with DFU development.

1. Introduction

Diabetic foot ulcers (DFUs) are among the most common complications affecting the lower extremities in diabetes mellitus patients [1]. These ulcers and complications thereof can cause high rates of morbidity and mortality among affected patients owing to associated angiopathy, oxidative microenvironmental damage, and repeated bacterial infections [2]. While substantial progress has been made in the treatment of DFUs in recent years, a large proportion nonetheless develops into chronic wounds through processes that are ultimately irreversible [3]. It is thus essential that the molecular mechanisms governing DFU development be clarified in order to aid in the prevention and treatment of these debilitating wounds [4].

Identifying genetic markers associated with DFU has the potential to guide the design of novel treatments while simultaneously elucidating the etiological basis for this condition [5]. Microarrays are commonly used to conduct large-scale bioinformatics studies aimed at simultaneously clarifying the relationship between multiple different genes and a given disease [6, 7]. One recent integrated bioinformatics analysis highlighted a role for the MAPK signaling pathway in DFU development [8]. Microarrays have also identified estrogen receptor 1, matrix metalloproteinase-2, and bone morphogenetic protein-4 as DFU-specific proteins [9-11]. In previous reports, differentially expressed genes (DEGs) associated with DFU progression have been attributed to a range of molecular functions, biological processes, and cellular structures. [12, 13]

In the present report, DFU-related DEGs were identified by analyzing previously published datasets containing DFU and normal tissue samples. We then conducted functional enrichment and protein-protein interaction (PPI) network analyses aimed at elucidating the mechanisms whereby these genes interact and cooperate to drive DFU development. Together, these results have the potential to clarify novel DFU-related biomarkers and to offer new insight regarding the molecular basis for this debilitating condition.

2. Materials and Methods

2.1. Dataset Selection

Microarray gene expression data of interest were downloaded from the Gene Expression Omnibus (GEO, http://www.ncbi.nlm.nih.gov/geo) database, which compiles a range of different high-throughput sequencing and microarray-based datasets. We searched this database for studies comparing DM2 and normal tissue samples and then downloaded the resultant files for analysis.

2.2. DEG Identification

Genes that were differentially expressed in DFU samples were identified using the default settings of the GEO2R tool (https://www.ncbi.nlm.nih.gov/geo/geo2r/), with DEGs being those genes with a P < 0.05 and a ∣logFC | >1. GEO2R was additionally used to construct volcano and Meandiff plots. Additionally, log2-transformed mRNA expression data were arranged into heatmaps using the “pheatmap” R package, while DEGs that were shared among datasets were determined using the Venn online tool (http://bioinformatics.psb.ugent.be/webtools/Venn/).

2.3. Functional Enrichment Analyses

The biological functions of identified DEGs of interest were assessed using the Database for Annotation, Visualization, and Integrated Discovery version (DAVID) Bioinformatics Resources (v6.8). Briefly, shared DEGs were imported into DAVID, and GO and KEGG enrichment analyses were then conducted. For GO analyses, enriched biological processes (BPs), molecular functions (MFs), and cellular components (CCs) were assessed. The “GOplot” R package was used to visualize the results of these enrichment analyses.

2.4. Protein-Protein Interaction (PPI) Network Analyses

To understand interactions among DEGs, PPI networks were constructed by importing up- and downregulated DEGs into the Search Tool for the Retrieval of Interacting Genes (STRING), with those interactions with a combined score >0.5 being used for network construction. Cytoscape (v 3.7.2) was used to visualize the network, while the cytoHubba plugin was used to rank genes within this network based upon their degree centrality values. Hub genes were considered to be those with the top 10 highest degree values.

3. Results and Discussion

3.1. DEG Identification

The GSE7014 and GSE29221 gene expression datasets were downloaded from the GEO database obtained from the GEO database. The GSE7014 dataset included 20 DM2 biopsy samples and 6 biopsy samples from normal individuals, whereas the GSE29221 dataset included 12 DM2 biopsies and 12 biopsies from normal individuals. In total, 1192 DEGs were identified in the GSE7014 dataset (900 upregulated, 292 downregulated), while 1177 were identified in the GSE29221 dataset (257 upregulated, 919 downregulated). The top 80 DEGs with the highest P values are presented in Figure 1. In total, 57 upregulated DEGs and 18 downregulated DEGs were shared between these two datasets, as identified through Venn diagram analyses (Figure 2).

Figure 1

Detection of differentially expressed genes (DEGs) in the GSE7014 and GSE29221 datasets. (a) An expression heat map of the top 80 DEGs in the GSE7014 dataset, as determined based upon P values. (b) A volcano plot corresponding to the GSE7014 dataset. (c) A Meandiff plot for the GSE7014 dataset. (d) An expression heat map of the top 80 DEGs in the GSE29221 dataset, as determined based upon P values. (e) A volcano plot corresponding to the GSE29221 dataset. (f) A Meandiff plot for the GSE29221 dataset.

Figure 2

Identification of shared DEGs. (a) DEGs upregulated in both the GSE7014 and GSE29221 datasets. (b) DEGs downregulated in both the GSE7014 and GSE29221 datasets.

3.2. Pathway Enrichment Analyses

GO analyses revealed these DEGs to be enriched in biological processes including sarcomere organization, muscle filament sliding, and the regulation of cardiac conduction, molecular functions including a structural constituent of muscle, protein binding, and calcium ion binding, and cellular components including Z disc, myosin filament, and M band. These DEGs were also enriched in KEGG pathways including the adrenergic signaling in cardiomyocytes, dilated cardiomyopathy, and tight junction pathways (Table 1 and Figure 3). Enrichment results pertaining to these analyses are compiled in Figure 4.

Table 1

Functional and pathway enrichment analyses for module genes. The top 3 terms were selected based upon p value rankings when >3 terms were enriched for a given category.

A, biological processes
Term	Name	Count	P value	Genes
GO:0045214	Sarcomere organization	7	1.3E-9	FHOD3, MYH3, ACTN2, CASQ2, CAPN3, LDB3, CASQ1
GO:0030049	Muscle filament sliding	4	4.4E-4	MYH3, ACTN2, MYL3, DMD
GO:1903779	Regulation of cardiac conduction	4	1.4E-3	PLN, CASQ2, ATP2B2, CASQ1
B, molecular functions
Term	Name	Count	P value	Genes
GO:0008307	Structural constituent of muscle	9	4.8E-12	MYOM1, PDLIM3, ACTN2, MYOT, MYL3, CAPN3, NEXN, DMD, MYOM2
GO:0005515	Protein binding	50	5.9E-4	FHOD3, MYOM1, BTG1, LGALSL, LDB3, SAT1, N4BP2L2, HK2,MYOM2, MED14, JPH2, XPO4, RASSF5, CAPN3, KIF1B, CTSC,ACTN2, IGFBP3, MYH1, PPP1R3C, MYOT, PRR16, CASQ2, TKT, FKBP3, CAMK2B, USP54, MGST1, THBS1, GTF2E2, PDLIM3, UGP2, PLN,CMYA5, KCNN2, DMD, PPARGC1A, S100A11, MPZL2, LGI1, MYH7B,DTNA, FAM46C, AGL, DDIT4L, ATP2B2, PPP2R3A, EFNA1, PKIA, CUTC
GO:0005509	Calcium ion binding	10	2.2E-3	ACTN2, MYL3, CASQ2, CAPN3, ATP2B2, PPP2R3A, CASQ1, THBS1, PLCD4, S100A11
C, cellular component
Term	Name	Count	P value	Genes
GO:0030018	Z disc	11	2.5E-11	FHOD3, JPH2, PDLIM3, ACTN2, MYOT, CASQ2, CAPN3, NEXN, LDB3, KCNN2, DMD
GO:0032982	Myosin filament	4	2.9E-5	MYH1, MYH7B, MYH3, MYOM2
GO:0031430	M band	4	9.1E-5	MYOM1, CMYA5, MYOM2, MYOM3
D, KEGG pathway
Term	Name	Count	P value	Genes
hsa04261	Adrenergic signaling in cardiomyocytes	5	4.5E-3	CAMK2B, PLN, MYL3, ATP2B2, PPP2R3A
Hsa05414	Dilated cardiomyopathy	4	8.0E-3	PLN, SGCD, MYL3, DMD
hsa04530	Tight junction	4	8.8E-3	MYH1, MYH7B, MYH3, ACTN2
KEGG: Kyoto Encyclopedia of Genes and Genomes.

Figure 3

GO term enrichment analysis results. (a) Z-score results for the top 6 GO terms, including the top 2 BPs, CCs, and MFs. (b) Enrichment results for DEGs and the top 6 GO terms. Z-scores were defined as follows: (upregulated genes–downregulated genes)/total genes.

Figure 4

KEGG pathway enrichment results. (a) Relationships between DEGs and the top 5 enriched KEGG pathways. (b) Cluster plots corresponding to DEGs and the top 5 enriched KEGG pathways.

3.3. PPI Network Construction and Hub Gene Identification

The STRING database was next used to construct a DEG PPI network (Figure 5), and the top 10 hub genes therein with the highest degree values were determined using Cytoscape v. 3.7.2. These hub genes were MYL3, ACTN2, DMD, PDLIM3, LDB3, MYH1, MYOM2, MYOT, CASQ2, and CAPN3 (Table 2 and Figure 6).

Figure 5

A DEG PPI network constructed using the STRING database.

Table 2

Degree of top 10 genes in top module.

Gene ID	Gene name	Degree
MYL3	Myosin Light Chain 3	31	Up
ACTN2	Actinin Alpha 2	30	Up
DMD	Dystrophin	26	Up
PDLIM3	PDZ And LIM Domain 3	24	Up
LDB3	LIM Domain Binding 3	24	Up
MYH1	Myosin Heavy Chain 1	22	Up
MYOM2	Myomesin 2	22	Up
MYOT	Myotilin	21	Up
CASQ2	Calsequestrin 2	21	Up
CAPN3	Calpain 3	21	Up

Up.

Figure 6

Hub gene identification. (a) A DEG PPI network constructed using Cytoscape, with upregulated and downregulated genes being shown in red and green, respectively. (b) The top 10 genes with the highest degree values were identified using CytoHubba. These genes were ranked in descending degree order from red to orange to yellow.

4. Discussion

Aberrant gene expression is closely linked to a range of pathological conditions, including DFU. Key driver genes linked to the onset and progression of this condition, however, remain to be fully clarified. In this study, we identified 900 upregulated and 292 downregulated DFU-related DEGs in the GSE7014 dataset, as well as 257 upregulated and 919 downregulated DFU-related DEGs in the GSE29221 dataset. These genes were associated with the adrenergic signaling in cardiomyocytes, dilated cardiomyopathy, and tight junction pathways. We were further able to identify 10 hub genes associated with DFU, including MYL3, ACTN2, DMD, PDLIM3, LDB3, MYH1, MYOM2, MYOT, CASQ2, and CAPN3.

DFU and other chronic wounds are associated with well-characterized morphological changes, but the underlying cellular and molecular biomarkers that drive these tissue changes remain poorly understood [14]. Changes in mRNA expression are valuable biomarkers that are well known to play a role in the development of diabetes-related diseases [15, 16]. For example, one prior study identified Prenylcysteine oxidase 1 (PCYOX1), beta-ala-his dipeptidase (CNDP1), and extracellular matrix protein 1 (ECM1) as valuable diagnostic biomarkers associated with the incidence of gestational diabetes [17]. Saik et al. further found the JAK-STAT, MAPK, and protein kinase B signaling pathway to be closely linked to diabetes complications and hypoxia responses [18]. The GO and KEGG analyses conducted in this study suggested the top DEGs to be enriched in sarcomere organization, muscle filament sliding, and the regulation of cardiac conduction, potentially playing a role in regulating angiogenesis. These genes were also enriched in molecular functions including a structural constituent of muscle, protein binding, and calcium ion binding, and in the adrenergic signaling in cardiomyocytes, dilated cardiomyopathy, and tight junction pathways, suggesting a potential role for the activation of inflammatory responses in DFU.

DFU-related hub genes identified in this study included MYL3, ACTN2, DMD, PDLIM3, LDB3, MYH1, MYOM2, MYOT, CASQ2, and CAPN3, all of which were involved in the top 5 KEGG pathways with the smallest P values. Degree centrality corresponds to the relationship between a given node and all other nodes in the network, while closeness centrality denotes the degree of closeness between a node and all other nodes in the network, and betweenness centrality measures the frequency with which a given node serves as the shortest bridge between two other nodes.

As most of the genes identified in this study have not previously been reported to be related to DFUs, there is a clear need to verify the functional importance and mechanistic roles of these genes in this pathological context. In addition, in vitro studies of human skin fibroblasts and human umbilical vein endothelial cells are warranted to explore the molecular mechanisms whereby these genes shape DFU development. The development of mice of other animal models in which these genes are conditionally knocked out may further aid in efforts to elucidate their functions as regulators of these debilitating chronic wounds.

5. Conclusions

In summary, the results of these bioinformatics analyses highlight novel mechanisms and important hub genes which may contribute to DFU development. However, further research is essential to better clarify the regulatory roles of these genes in order to firmly establish their value as clinical biomarkers and/or therapeutic targets.