Abstract

Purpose

   To investigate the molecular mechanism and search for candidate
   biomarkers in the gene expression profile of patients with diabetic
   peripheral neuropathy (DPN).

Methods

   Differentially expressed genes (DEGs) of progressive vs non-progressive
   DPN patients in dataset [28]GSE24290 were screened. Functional
   enrichment analysis was conducted, and hub genes were extracted from
   the protein–protein interaction network. The expression level of hub
   genes in serum samples in another dataset [29]GSE95849 was obtained,
   followed by the ROC curve analysis.

Results

   A total of 352 DEGs were obtained from dataset [30]GSE24290. They were
   involved in 14 gene ontology terms and 10 Kyoto Encyclopedia of Genes
   and Genomes pathways, mainly related to lipid metabolism. Eight hub
   genes (LEP, APOE, ADIPOQ, FABP4, CD36, GPAM, CIDEC, and PNPLA4) were
   revealed, and their expression level was obtained in dataset
   [31]GSE95849. The receiver operating characteristic curve analysis
   indicated that CIDEC (AUC=1), APOE (AUC=0.833), CD36 (AUC=0.803), and
   PNPLA4 (AUC=0.861) might be candidate serum biomarkers of DPN.

Conclusion

   Lipid metabolism of Schwann cells might be inhibited in progressive
   DPN. CIDEC, APOE, CD36, and PNPLA4 might be potential predictive
   biomarkers in the early DPN diagnosis of patients with DM.

   Keywords: differentially expressed genes, functional enrichment
   analysis, demyelination, lipid metabolism, 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.[32]^1 Approximately 50% of those with diabetes may develop a
   diabetic peripheral neuropathy (DPN). The number will only increase as
   the diabetes epidemic grows.[33]^2 DPN is characterized by pain,
   paraesthesia, and sensory loss.[34]^3 Patients complaint unbearable
   lancinating, tingling and burning sensation, even depression, anxiety,
   and sleep deprivation.[35]^4^,[36]^5 Moreover, insensitivity to trauma
   often results in foot ulcers which can lead to some levels of
   amputation.[37]^6 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.[38]^7 Various
   symptomatic medications are nonspecific with serious side
   effects.[39]^8^,[40]^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.[41]^10 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.[42]^11 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 [43]GSE24290 deposited by Hur et al.[44]^11 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.[45]^12^,[46]^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.[47]^14^–[48]^16
   Then, pathview was used to describe the most important
   pathway.[49]^17^,[50]^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).[51]^19 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.[52]^20
   Connectivity degree analysis was performed and the most highly
   connected cluster was extracted from the PPI network through MCODE
   analysis.[53]^21

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 [54]GSE95849, which contains serum samples
   from six DM patients and six DPN patients.[55]^22 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).[56]^23
   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 [57]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 ([58]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 ([59]Figure 1A and [60]B and [61]Figure S2)

Figure S1.

   [62]Figure S1
   [63]Open in a new tab

   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 1.

   [64]Figure 1
   [65]Open in a new tab

   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 S2.

   Figure S2
   [66]Open in a new tab

   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
   ([67]Table 1). Kyoto Encyclopedia of Genes and Genomes (KEGG)
   enrichment analysis unveiled ten significantly enriched KEGG pathways
   with P-value <0.05 ([68]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-γ ([69]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 [70]Table 2.

Table 1.

   GO terms enrichment results of DEGs
   Term Count P-value Enriched genes
   Enriched from up-regulated genes
    GO:0007154~cell communication 3 0.0169 GJB6, GJB2, GJA3
    GO:0050896~response to stimulus 3 0.0253 CLRN1, KERA, RPGRIP1
   Enriched from down-regulated genes
    GO:0055089~fatty acid homeostasis 3 0.0031 APOE, MLXIPL, GPAM
    GO:0042593~glucose homeostasis 5 0.0042 LEP, CAV3, GPR21, MLXIPL,
   ADIPOQ
    GO:0006936~muscle contraction 5 0.0052 CAV3, MYH1, ACTA1, MYH13, MYH7
    GO:0042632~cholesterol homeostasis 4 0.0083 CAV3, MALRD1, APOE, FABP4
    GO:0019433~triglyceride catabolic process 3 0.0114 APOE, FABP4, PNPLA3
    GO:0051897~positive regulation of protein kinase B signaling 4 0.0173
   LEP, TNFAIP8L3, STOX1, IL26
    GO:0002027~regulation of heart rate 3 0.0193 CAV3, RYR2, MYH7
    GO:0043407~negative regulation of MAP kinase activity 3 0.0223 CAV3,
   APOE, ADIPOQ
    GO:0006810~transport 7 0.0259 CLCA2, AQP8, GRIK4, FABP4, CLVS1, CNGA3,
   TRPM1
    GO:0035338~long-chain fatty-acyl-CoA biosynthetic process 3 0.0304
   SCD, ACSBG2, ACOT6
    GO:0006635~fatty acid beta-oxidation 3 0.0331 LEP, ABCD2, ADIPOQ
    GO:0034220~ion transmembrane transport 5 0.0480 CLCA2, AQP8, GRIK4,
   ANO3, RYR2
   [71]Open in a new tab

   Notes: Significantly enriched GO terms with P-value <0.05 count ≥3 and
   were screened out. Five terms were directly associated with lipid
   metabolism.

   Abbreviations: GO, Gene Ontology; DEGs, differentially expressed genes.

Figure S3.

   [72]Figure S3
   [73]Open in a new tab

   The PPAR signaling pathway and DEGs.

   Abbreviation: DEGs, differentially expressed genes.

Table 2.

   Information of KEGG analysis of dataset [74]GSE24290
   Term Count P-value Enriched genes Fold enrichment
   hsa03320:PPAR signaling pathway 7 <0.0001 CD36, OLR1, SCD, ACSBG2,
   FABP4, ADIPOQ, MMP1 11.0569
   hsa04923:Regulation of lipolysis in adipocytes 5 0.0002 PTGER3, NPY,
   PDE3B, FABP4, PRKG2 9.4492
   hsa00982:Drug metabolism – cytochrome P450 5 0.0004 CYP2A6, ADH1A,
   GSTM5, UGT2B28, MGST1 7.7817
   hsa04920:Adipocytokine signaling pathway 5 0.0004 LEP, CD36, NPY,
   ACSBG2, ADIPOQ 7.5593
   hsa00980:Metabolism of xenobiotics by cytochrome P450 5 0.0005 CYP2A6,
   ADH1A, GSTM5, UGT2B28, MGST1 7.1507
   hsa05204:Chemical carcinogenesis 5 0.0006 CYP2A6, ADH1A, GSTM5,
   UGT2B28, MGST1 6.6144
   hsa04060:Cytokine-cytokine receptor interaction 8 0.0007 LEP, IL20RB,
   INHBE, GDF5, IL26, EDAR, CCL27, CCL26 3.4841
   hsa00830:Retinol metabolism 4 0.0213 CYP2A6, SDR16C5, ADH1A, UGT2B28
   6.6144
   hsa04152:AMPK signaling pathway 5 0.0271 LEP, CD36, SCD, GYS2, ADIPOQ
   4.3021
   hsa04024:cAMP signaling pathway 6 0.0364 PTGER3, NPY, PDE3B, RYR2,
   CNGA3, HTR1D 3.2070
   [75]Open in a new tab

   Notes: KEGG biological pathway enrichment analysis found that the PPAR
   signaling pathway (P-value <0.0001, count =7, and fold enrichment
   =11.0569) was the most important one among the ten significantly
   enriched pathways.

   Abbreviation: KEGG, Kyoto Encyclopedia of Genes and Genomes.

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 ([76]Figure 2A).
   Subsequently, we analyzed the PPI network and extracted the most highly
   connected cluster by MCODE plug-in in Cytoscape ([77]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.

   [78]Figure 2
   [79]Open in a new tab

   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
   ([80]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.

   [81]Figure 3
   [82]Open in a new tab

   The ROC curves of hub genes in [83]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^2,
   while samples in the non-progressor group lost ≤100 fibers/mm^2 over 52
   weeks.[84]^11

   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” ([85]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.[86]^24 After being damaged by hyperglycemia, the balance
   of Schwann cells de-differentiation and re-differentiation can be
   destroyed.[87]^25 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.[88]^26 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.[89]^27 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.[90]^28^–[91]^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.[92]^31^–[93]^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.[94]^36^–[95]^39 PPAR-γ agonist may lead to upregulation of
   CIDEC thereby increasing lipid accumulation.[96]^40^,[97]^41 Although
   the expression of CIDEC proven to positively correlate with the
   development of insulin sensitivity in obese people,[98]^42 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.[99]^43

   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 [100]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 [101]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 [102]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 ([103]GSE24290 and [104]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.

Acknowledgments