Gene expression profiling reveals candidate biomarkers and probable molecular mechanism in diabetic peripheral neuropathy

Introduction

It was estimated that 415 million people aged 20–79 years suffered from diabetes in 2015, and the number was predicted to rise to 642 million by 2040.¹ Approximately 50% of those with diabetes may develop a diabetic peripheral neuropathy (DPN). The number will only increase as the diabetes epidemic grows.² DPN is characterized by pain, paraesthesia, and sensory loss.³ Patients complaint unbearable lancinating, tingling and burning sensation, even depression, anxiety, and sleep deprivation.^4,5 Moreover, insensitivity to trauma often results in foot ulcers which can lead to some levels of amputation.⁶ Several attempts have been made to treat this disease; however, none of the pharmacotherapy have proven to be effective in altering the progressive course to date.⁷ Various symptomatic medications are nonspecific with serious side effects.^8,9 The treatment dilemma reflects our present knowledge for the pathogenesis of DPN is far from to be clear.

Demyelination is known as an early pathological feature in DPN, and it precedes the degeneration of axon.¹⁰ Therefore, studying the molecular mechanism associated with demyelination may conduce to a better understanding of DPN and finding the predictive biomarkers of this disease. In 2010, Hur et al, analyzed genes that differed between DPN progressors and non-progressors, which is classified by the loss of myelinated fiber density of sural nerves.¹¹ Based on the microarray data, the present study analyzed the differentially expressed genes (DEGs) via function enrichment analysis and topological approaches. We found several hub genes according to the most significant cluster of protein–protein interaction (PPI) network. The expression level of each hub gene was obtained in serum samples of another dataset, followed by receiver operating characteristic (ROC) curve study. The results may reveal the molecular mechanism of demyelination in DPN and provide potential serum biomarkers.

Materials and methods

Data preprocessing and DEGs screening

From the Gene Expression Omnibus database, we collected the gene expression data of GSE24290 deposited by Hur et al.¹¹ There were 35 specimens including progressors and non-progressors sural nerve biopsies for analysis. The raw expression data underwent background correction, normalization, and summarization using the robust multi-array average (RMA) algorithm in oligo.^12,13 The limma package in R was applied to identify DEGs between two groups. In this analysis, P-value <0.05 and log|FC| >0.5 were used as the cutoff criteria.

Functional enrichment and PPI network analysis

To investigate the main functional pathways of DPN, we submitted the DEGs to Database for Annotation, Visualization, and Integrated Discovery (DAVID), which was used to perform the Gene Ontology (GO) analysis and KEGG pathways enrichment analysis of DEGs.^14–16 Then, pathview was used to describe the most important pathway.^17,18 Criteria for this step were set as P-value<0.05 and gene counts ≥3. The functional protein interactions of DEGs and the encoding proteins were predicted using the Search Tool for the Retrieval of Interacting Genes (STRING).¹⁹ Cytoscape 3. 6. 1 was used for visualization and to calculate the properties of the PPI network. Subsequently, the Network Analyzer plug-in of Cytoscape was utilized to analyze the topology properties of the network.²⁰ Connectivity degree analysis was performed and the most highly connected cluster was extracted from the PPI network through MCODE analysis.²¹

ROC curve analysis

The genes that constituted the most highly connected cluster in PPI network were considered to be hub genes. Then, we downloaded the gene expression profiling data of GSE95849, which contains serum samples from six DM patients and six DPN patients.²² We obtained the expression level of hub genes and used the pROC package in R software to draw the ROC curves and calculate the area under curve (AUC).²³ Larger AUC value means the gene can well distinguish DPN from the DM patient samples. The diagnosis effect of hub genes was further investigated according to the AUC value.

Results

DEGs identified between progressive and non-progressive DPN patients

The raw gene expression data of GSE24290 were normalized by RMA. The boxplots of the intensity of all samples demonstrated that the expression values of each sample were close to the same after normalization (Figure S1). A total of 352 DEGs between non-progressive and progressive DPN patients were obtained with the criteria of P-value <0.05 and log|FC|>0.5. The 142 up-regulated genes and the 210 down-regulated genes were shown in the heat map and the volcano plot (Figure 1A and B and Figure S2)

Figure 1 The heat map and the volcano plot of differentially expressed genes. Notes: (A) The heat map of the differentially expressed genes in progressive DPN group vs non-progressive DPN group. The horizontal axis represents the 35 samples, and the vertical represents the top50 up-regulated and top50 down-regulated genes sorted by P-value. Up-regulated genes are shown by warm colors. Down-regulated genes are shown by cool colors. (B) The volcano plot of the differentially expressed genes. Red dots on the right indicate up-regulated genes, and blue dots on the left indicate downregulation. Gray dots indicate genes with no statistically significant difference. (C) The bubble map of KEGG pathway analysis. The horizontal axis represents the fold enrichment of pathways, and the vertical represents pathway names. Size of bubbles represents the number of genes, and the shade of color depends on the P-value.Abbreviations: DPN, diabetic peripheral neuropathy; DEG, differentially expressed genes; FC, fold-change; KEGG, Kyoto Encyclopedia of Genes and Genomes.

Figure S1 The boxplots of sample data before and after normalization. Notes: The lateral axis represented 35 samples from 18 progressive and 17 non-progressive DPN patients. The longitudinal axis represents expression levels. The horizontal line in the middle of post represents the expression levels of each sample. Each sample was close to the same following normalization.Abbreviation: DPN, diabetic peripheral neuropathy.

Figure S2 The heat map of all 352 DEGs.Abbreviation: DEGs, differentially expressed genes.

Functional annotation and enrichment analysis of the DEGs

To reveal further insights into the biological functions of DEGs, functional enrichment analyses were performed using DAVID. For the GO analysis, we focused on the categories of biological processes (BP) and set the criteria with P-value <0.05 as significantly enriched GO terms. As a result, DAVID identified two terms significantly enriched from 142 up-regulated genes and 12 terms significantly enriched from 210 down-regulated genes. These terms demonstrated that the gene expression differences of nerve samples in progressive DPN were associated with “fatty acid homeostasis”, “glucose homeostasis”, and other BP (Table 1). Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis unveiled ten significantly enriched KEGG pathways with P-value <0.05 (Figure 1C). Among these pathways, the PPAR signaling pathway was the most important one according to the P-value, gene counts, and fold enrichment. The DEGs in the PPAR pathway were concentrated in PPAR-γ (Figure S3). There were seven DEGs genes that participated in this pathway (CD36, OLR1, SCD, ACSBG2, FABP4, ADIPOQ, MMP1). We listed other useful information of these KEGG pathways in Table 2.

Table 1 GO terms enrichment results of DEGs

Table 2 Information of KEGG analysis of dataset GSE24290

Figure S3 The PPAR signaling pathway and DEGs.Abbreviation: DEGs, differentially expressed genes.

PPI network clusters analysis and selection of hub genes

PPI network was constructed involving 73 nodes (DEGs) and 132 edges. The nodes represented the proteins expressed by DEGs and the edges between two nodes means the physical interactions. The connectivity degree is an important parameter and the high connectivity degree indicated the protein interacted with more surrounding proteins and play a more important role. The nodes with higher connectivity degree were shown as larger sizes with red or orange color (Figure 2A). Subsequently, we analyzed the PPI network and extracted the most highly connected cluster by MCODE plug-in in Cytoscape (Figure 2B). Genes in this cluster were at the core of the whole network, including LEP, APOE, ADIPOQ, FABP4, CD36, PNPLA4, GPAM, and CIDEC. Hence, we consider the eight genes as the hub genes for further analysis, and they were all down-regulated as demyelination.

Figure 2 The protein–protein interaction network and the most highly connected cluster. Notes: (A) The protein–protein interaction network consists of 73 nodes and 132 edges. Color and size represent the connectivity degree of nodes. (B) The most highly connected cluster is composed of eight hub genes: LEP, APOE, ADIPOQ, CD36, FABP4, CIDEC, GPAM, and PNPLA4.

ROC curves and candidate biomarkers

The ROC curve analysis was performed by the pROC package in R (Figure 3). The AUCs of each gene were more than 0.500. Among these hub genes, CIDEC was the outstanding one with the AUC =1.000, which represented it might have great value for the diagnosis of DPN in patients with DM. Besides, APOE (AUC=0.833), CD36 (AUC=0.803), and PNPLA4 (AUC=0.861) could also be predictive biomarkers. The AUCs of other hub genes were less than 0.8000.

Figure 3 The ROC curves of hub genes in GSE95849. Notes: CIDEC, PNPLA4, APOE, and CD36 were four genes with the AUC>0.8000.Abbreviations: ROC, receiver operating characteristic; AUC, the area under curve.

Discussion

The mechanism producing DPN is multifactorial and extremely complex. To further understand the molecular mechanism and search for novel serum biomarkers for DPN, we analyzed two different datasets of expression profile by bioinformatic approaches in the current study. We first compared the microarray data of two groups of sural nerve samples from DPN progressors and non-progressors. They were divided by the level of demyelination. Samples in the progressor group lost ≥500 fibers/mm², while samples in the non-progressor group lost ≤100 fibers/mm² over 52 weeks.¹¹

Totally 352 DEGs were identified for following functional enrichment analysis and PPI network analysis. The results demonstrated the function of DEGs was closely related to lipid metabolism. In GO analysis, most of the terms were derived from down-regulated DEGs. We only focused on the most prominent parts of these terms. It was not difficult to find that five terms were directly associated with lipid metabolism in total 14 GO terms of BP. They were “fatty acid homeostasis”, “cholesterol homeostasis”, “Adipocytokine signaling pathway”, “triglyceride catabolic process”, “long-chain fatty-acyl-CoA biosynthetic process”, and “fatty acid beta-oxidation” (Table 1). The changes were likely to occur in Schwann cells. In the peripheral nerve system, myelin results from the circumferential wrapping of the Schwann cell plasma membrane, which is enriched in lipid-like glycosphingolipids, saturated long-chain fatty acids and, particularly, cholesterol.²⁴ After being damaged by hyperglycemia, the balance of Schwann cells de-differentiation and re-differentiation can be destroyed.²⁵ As the fatty acid homeostasis, cholesterol homeostasis was changed and the triglyceride catabolic process, long-chain fatty-acyl-CoA biosynthetic process was inhibited, re-differentiated mature Schwann cells might be difficult to generate new myelin, which might be an important cause of DPN progression.

The KEGG analysis further revealed changes in the signaling pathways in progressive DPN. Three important signaling pathways were associated with lipid metabolism. They were “PPAR signaling pathway”, “regulation of lipolysis in adipocytes”, and “adipocytokine signaling pathway”. The “PPAR signaling pathway” is highlighted with smallest P-value and largest fold enrichment among them. Almost all DEGs involved in the PPAR signaling pathway were down-regulated, which was consistent with the study of Kim et al.²⁶ They demonstrated chronic high glucose inhibited the function of PPAR-γ due to the reduction of PPAR-γ binding to target genes in Schwann cells. Montani et al, described the endogenous fatty acid synthesis, which was potentially critical process of myelination, could trigger activation of the PPAR-γ transcriptional program in Schwann cells and the PPAR-γ agonist could partially rescue Schwann cell myelination in the setting of deficient endogenous fatty acid synthesis.²⁷ Our study suggested the inhibition of PPAR-γ signaling pathway caused by hyperglycemia might be crucial in the progression of DPN.

The PPI network contained multiple clusters. We applied the MCODE in Cytoscape to extract the most significant cluster of PPI network. All the genes that make up this cluster play important roles in lipid metabolism. For example, CD36 facilitates cell membrane free fatty acid transport in adipocytes and ADIPOQ and LEP were secreted by adipocytes.^28–30 Accordingly, we assumed that Schwann cells could secrete some adipokines like adiponectin and take in lipid and free fatty acid via CD36. Uptake of palmitic acids, which is the most abundant plasma free fatty acid involved in inducing insulin resistance, proven to cause Schwann cells dysfunction and death.^31–35 Downregulation of CD36 in progressive DPN might reduce lipid uptake and influence myelination, but might play a protective role.

Furthermore, we suggested that CIDEC, PNPLA4, APOE, and CD36 might be used as potential molecules for liquid biopsy of DPN based on the results of ROC curve analysis. In particular, CIDEC might be an important molecule that has not been thoroughly studied but is of great value to the diagnosis of DPN. CIDEC protein is a member of the cell death-inducing DNA fragmentation factor-like effector family, which are crucial for multiple lipid metabolic pathways and lipid homeostasis.^36–39 PPAR-γ agonist may lead to upregulation of CIDEC thereby increasing lipid accumulation.^40,41 Although the expression of CIDEC proven to positively correlate with the development of insulin sensitivity in obese people,⁴² it has not been studied within the context of DPN by now. The only investigation of CIDEC associated with diabetic complications is that silence of CIDEC may partially reverse diabetic pulmonary vascular in type 2 diabetes.⁴³

The present study re-analyzed the two publicly available microarray gene expression profiling via different methods and arrived at different conclusions. We constructed the PPI network of DEGs of the dataset GSE24290 and extracted the most important cluster via MCODE. This topological analysis differs from gene co-citation analysis in primary publication. Combining with our results from analysis of the dataset GSE95849, we proposed our new point that lipid metabolism of Schwann cells might be inhibited in progressive DPN and the inhibition of PPAR-γ signaling pathway might be crucial in the pathogenesis of the disease. In addition, in the primary publication of the dataset GSE95849, the research found the downregulation of the neurotrophin-MAPK signaling pathway may be crucial for DPN pathogenesis, while our analysis only focused on the expression level of several selected genes in serum samples of the dataset.

Conclusion

The present study aimed to investigate the molecular mechanism in gene expression profiling in the pathogenesis of DPN. Totally 352 DEGs and eight hub genes were screened via this bioinformatic approaches of two microarray datasets (GSE24290 and GSE95849). Our salient findings were that lipid metabolism of Schwann cells might be inhibited in progressive DPN and CIDEC, APOE, CD36, PNPLA4 were identified as candidate predictive biomarkers in the early DPN diagnosis of patients with DM. However, there were some limitations in the present study such as small sample size and lack of verification test. Further basic experiments with large sample size are needed to validate our results.