Skip to main content

Genome-wide profiling of DNA methylation and gene expression identifies candidate genes for human diabetic neuropathy

  • The Correction to this article has been published in Clinical Epigenetics 2020 12:130

Abstract

Background

Diabetic peripheral neuropathy (DPN) is the most common complication of type 2 diabetes (T2D). Although the cellular and molecular mechanisms of DPN are poorly understood, we and others have shown that altered gene expression and DNA methylation are implicated in disease pathogenesis. However, how DNA methylation might functionally impact gene expression and contribute to nerve damage remains unclear. Here, we analyzed genome-wide transcriptomic and methylomic profiles of sural nerves from T2D patients with DPN.

Results

Unbiased clustering of transcriptomics data separated samples into groups, which correlated with HbA1c levels. Accordingly, we found 998 differentially expressed genes (DEGs) and 929 differentially methylated genes (DMGs) between the groups with the highest and lowest HbA1c levels. Functional enrichment analysis revealed that DEGs and DMGs were enriched for pathways known to play a role in DPN, including those related to the immune system, extracellular matrix (ECM), and axon guidance. To understand the interaction between the transcriptome and methylome in DPN, we performed an integrated analysis of the overlapping genes between DEGs and DMGs. Integrated functional and network analysis identified genes and pathways modulating functions such as immune response, ECM regulation, and PI3K-Akt signaling.

Conclusion

These results suggest for the first time that DNA methylation is a mechanism regulating gene expression in DPN. Overall, DPN patients with high HbA1c have distinct alterations in sural nerve DNA methylome and transcriptome, suggesting that optimal glycemic control in DPN patients is an important factor in maintaining epigenetic homeostasis and nerve function.

Background

Type 2 diabetes (T2D) accounts for 90 to 95% of the 463 million people suffering from diabetes worldwide [1] and is associated with macro- and microvascular complications. Of these, microvascular complications are more common and affect multiple tissues including the eye, kidney, and nerve [2, 3]. Diabetic neuropathy affects the nerve and is the most prevalent microvascular complication that can present in multiple forms, the most common being diabetic peripheral neuropathy (DPN) [4]. DPN affects up to 60% of T2D patients and is characterized by distal-to-proximal nerve damage that results in reduced sensation, chronic pain, increased infection risk, and recurrent foot ulceration that can lead to lower limb amputations [4,5,6]. Despite the enormity of the problem, current treatments often fail to slow or reverse DPN progression in T2D [7]. Therefore, understanding the pathogenesis of DPN is critical for developing mechanism-based therapies.

The risk of developing T2D and DPN is determined by both genetic and lifestyle factors [3, 8]. Transcriptome profiling has provided insight into disease pathogenesis by identifying numerous genes and pathways implicated in DPN, including inflammation [9,10,11,12], oxidative stress [10, 12, 13], lipid metabolism [14, 15], and mitochondrial dysfunction [11, 16]. However, many of these previous analyses were either performed in animal models or using microarrays, which have several limitations, including a narrow dynamic range, low sensitivity, and dependency on predefined transcripts [17,18,19]. Next-generation RNA-sequencing (RNA-seq) is considered a potential alternative to microarrays because it is unbiased and more sensitive [18, 20]. However, RNA-seq has not yet been used to characterize the transcriptomic changes in human DPN.

In addition to genetic factors, epigenetic mechanisms, such as DNA methylation, histone modifications, chromatin remodeling, and RNA editing and biogenesis have recently emerged as a potential link between gene expression and environmental factors [21]. DNA methylation refers to the reversible attachment of a methyl group to a cytosine within cytosine–phosphate–guanine (CpG) dinucleotides [22]. In differentiated cells, DNA methylation contributes to the maintenance of normal DNA structure, chromosome stability, and gene regulation [23]. DNA methylation regulates gene expression without altering the underlying DNA sequence and is of particular interest because of its emerging role in T2D and its complications [24,25,26,27]. We recently showed that aberrant DNA methylation is involved in nerve degeneration in T2D and DPN in a small cohort of patients [24]. Specifically, our results highlighted the role of DNA methylation in regulating pathways previously shown to be implicated in DPN pathogenesis, including axon guidance, glycerophospholipid metabolism, and MAPK signaling. However, much less is known about the impact of differential DNA methylation on gene expression in DPN and how the interaction between genetic and epigenetic mechanisms may affect biological pathways during DPN pathogenesis.

In this study, we expanded on our previous findings by conducting a comprehensive systems biology analysis that combines global DNA methylome and genome-wide transcriptome profiling in human sural nerves of T2D patients with DPN. By integrating epigenomic and transcriptomic data, we found associations between hemoglobin A1c (HbA1c), DNA methylation, and differential gene expression patterns within pathways known to play a key role in DPN pathogenesis. These findings support the involvement of epigenetic dysregulation in DPN and offer candidate targets for developing therapeutic strategies for DPN.

Results

Study population

To identify specific mechanisms that are differentially transcribed in T2D and DPN, we determined gene expression profiles using RNA-seq in sural nerve biopsies obtained from DPN patients from a double-blind placebo-controlled clinical trial of a candidate treatment that proved ineffective [28, 29]. Sural nerve biopsies (n = 78) from T2D patients with DPN were initially classified based on changes in myelinated fiber density (MFD) over 52 weeks [30]. Genome-wide DNA methylation profiling was also performed on the same cohort, using reduced representation bisulfite sequencing (RRBS). Patient clinical characteristics are available in Additional file 2: Table S1, including HbA1c, triglyceride and cholesterol levels, and neuropathy phenotyping by MFD counts and O’Brien neuropathy score.

RNA-seq and RRBS quality assessment and alignment summary

All samples had high-quality RNA-seq reads with an average of 96% reads per sample that survived trimming with Trimmomatic (Additional file 2: Table S2). Approximately 77% of the total reads uniquely aligned to the hg19 reference human genome, except one sample (ID 63096) with a unique mapping rate of 33%, which was excluded from subsequent analyses. In total, an average of 90 million paired-end and single-end reads per sample was generated by RNA-seq. To identify genes and pathways under epigenetic control in human DPN, we next examined genome-wide DNA methylation using RRBS profiling. Briefly, over 98% were high-quality reads, and approximately 52 million reads per sample could be uniquely mapped (55%~68%) to hg19 (Additional file 2: Table S2). CpG sites on X and Y chromosomes were excluded, resulting in a total of 1,602,934 CpG sites identified for differential methylation analysis. In total, an average of 45 million pairs of 50-bp paired-end reads and 90 million 50-bp single-end reads were generated for RNA-seq and RRBS, respectively.

Data-driven clustering analysis

We first examined the overall gene expression similarity among the samples using hierarchical clustering analysis. Cluster dendrogram (Additional file 1: Figure S1A) defined three groups of samples, which did not correlate with the previous classification based on changes in MFD [30]. We then identified the clinical factors that were significantly associated with these transcriptomics data-driven groups by applying principal component analysis (Additional file 1: Figure S1B) and multifactorial logistic regression analysis. HbA1C levels were the only clinical factor that significantly differed across these groups (p = 0.04; Fig. 1). For downstream analyses, we focused on clusters with the highest (group 1, n = 21) and lowest (group 2, n = 32) HbA1c. Clinical characteristics and neuropathy measures for groups 1 and 2 are summarized in Table 1.

Fig. 1
figure1

Violin plot of % hemoglobin A1c (HbA1c) level distribution in groups 1 and 2. Each dot corresponds to the HbA1c (%) of a patient and color corresponds to the groups

Table 1 Clinical characteristics and neuropathy measures of groups 1 and 2

Differentially expressed genes (DEGs) analysis

A total of 998 genes were differentially expressed in group 1 (high HbA1c) relative to group 2 (low HbA1c) with a Benjamini-Hochberg adjusted p value < 0.01. Of the 998 DEGs, we found that 542 genes were significantly upregulated and 456 genes were significantly downregulated (Additional file 2: Table S3). The top 30 most upregulated and downregulated DEGs along with their annotation are presented in Table 2. Additionally, the top 50 most significant DEGs are labeled in an MA plot (Additional file 1: Figure S2) and gene expression levels of the 100 most significant DEGs in group 1 versus group 2 are displayed in a heatmap (Additional file 1: Figure S3). Interestingly, significantly upregulated genes in group 1 compared to group 2 included the long noncoding RNA NEAT1, which has been recently implicated in T2D and neurodegeneration [31, 32]. The immune response has also been repeatedly shown to be involved in T2D and neurodegeneration [9, 14], and downregulated DEGs included immune response genes CCDC80 and CD14.

Table 2 The 15 most upregulated and downregulated DEGs

Functional enrichment analyses of DEGs

To better understand biological functions and molecular pathways associated with the DEGs, we performed functional enrichment analysis using Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway, Gene Ontology (GO), Reactome, and Disease Ontology (DO) terms with an adjusted p value < 0.05 as the significance cutoff. Using KEGG pathway analysis, we found significant enrichment of pathways related to phagosomes as well as antigen processing and presentation (Fig. 2a, Additional file 2: Table S4). “Cell adhesion molecules (CAMs),” “Type I diabetes mellitus,” and “Neurotrophin signaling pathway” were also significantly enriched (Fig. 2a, Additional file 2: Table S4). Moreover, 688 GO terms were significantly overrepresented in group 1 relative to group 2. Further clustering based on gene content using Database for Annotation, Visualization and Integrated Discovery (DAVID) clustering [33] identified 24 significant GO clusters. Those included biological processes such as “cellular component organization,” “extracellular matrix organization,” and “regulation of immune response” (Fig. 2b, Additional file 2: Table S5). Within these pathways, top DEGs included members of the transforming growth factor-beta (TGF-β) family, previously associated with the development of DPN [34]. Additionally, pathways involving calcium signaling and homeostasis were dysregulated in group 1 compared to group 2, with significant enrichment in terms like “calcium-mediated signaling” and “calcium ion transport into cytosol” (Additional file 2: Table S5).

Fig. 2
figure2

Functional enrichment analysis of DEGs by KEGG and GO. The 20 most significantly enriched biological functions using KEGG (a) and GO (b) are illustrated in dot plots. Rich factor refers to the proportion of DEGs belonging to a specific term. Node size (gene number) refers to the number of DEGs within each term and node color indicates the level of significance (−log10  p value)

We also performed enrichment analysis using Reactome and DO, and found that 58 Reactome and 39 DO pathways were significantly altered (adjusted p value < 0.05, Additional file 1: Figure S4A-B, Additional file 2: Tables S6-7). Many of these pathways were similar to those identified by KEGG and GO analysis, including “adaptive immune system,” “extracellular matrix organization,” “axon guidance,” and “demyelinating disease.” Interestingly, DEGs were highly involved in neuronal and glial function as well as myelination according to DO, such as insulin-like growth factor I (IGF-I) (Additional file 2: Tables S6-7) and the antioxidant superoxide dismutase 2 (SOD2) (Additional file 2: Table S7) [35, 36].

DNA methylation changes

We found a total of 2066 differentially methylated CpG sites (DMCs) between group 1 and group 2, of which 1169 were hypomethylated and 897 were hypermethylated (Additional file 2: Table S8). The overall percentage of DMCs ranged between 45 and 50% (Fig. 3a), with no significant differences in methylation patterns across the chromosomes. We also determined the genomic distribution of altered DNA methylation and the distribution of DMCs in relation to CpG islands. Approximately 7% of DMCs were located in promoter regions, 37% in introns, 12% in exons, and ~ 44% in intergenic regions (Fig. 3b). Additionally, the majority of DMCs were in “non-CpG-rich regions” (denoted as “other” in the figure), while 34% of DMCs were located in CpG islands and shores (Fig. 3c).

Fig. 3
figure3

Distribution of differentially methylated CpGs (DMCs). a DMCs across chromosomes are depicted in a circular plot. Hyper- and hypomethylated CpGs are colored in red and blue, respectively, relative to group 2. The distributions of DMCs summarized based on genomic location (b) and relative to CpG islands (CpGi) (c)

To explore how DMCs induce functional changes in the genome, we mapped them to 1519 unique NCBI Reference Sequence Database (RefSeq) IDs (Additional file 2: Table S8). Only the genes with valid official gene symbols, referred to as differential methylated genes (DMGs) (n = 1489), were included in downstream analyses. A total of 929 DMGs (mapped from 1167 DMCs) exhibited hypomethylation, while 686 DMGs (mapped from 894 DMCs) exhibited hypermethylation in group 1 compared to group 2. Approximately 8% of DMGs (n = 126) were mapped from both hypermethylated and hypomethylated DMCs. With respect to the distance of DMCs from transcript start sites (TSSs), 201 DMCs corresponding to 153 DMGs were within a 5 kb distance upstream from their respective TSSs. The top 20 DMCs within the promoter regions (< 5 kb distance from TSS), along with their annotated genes (DMGs), are listed in Table 3 and include several microRNAs and coding genes.

Table 3 Top 20 DMCs within the promoter region

Functional enrichment analysis of DMGs

We next annotated DMGs for biological and functional enrichment using KEGG, GO, Reactome, and DO terms (Fig. 4, Additional file 1: Figure S5). A total of 24 KEGG pathways were significantly enriched and included “extracellular matrix (ECM)-receptor interaction,” “MAPK signaling pathway,” “ErbB signaling pathway,” and “axon guidance” (Fig. 4a, Additional file 2: Table S9). Significant DMGs within these pathways included structural genes and genes related to nerve damage, inflammation, and diabetes, such as ERBB4 and IL1B. GO analysis showed that 1519 DMGs were highly enriched in 929 biological processes, which were further clustered using DAVID clustering into 15 GO clusters with kappa > 0.5. Over-represented clusters included processes related to ECM remodeling such as cell adhesion as well as nervous system development (Fig. 4b, Additional file 2: Table S10), with significant DMGs such as SOX9 and Nr2f2 (Additional file 2: Table S10).

Fig. 4
figure4

Functional enrichment analysis of DMGs by KEGG and GO. The 20 most significantly enriched biological functions using KEGG (a) and GO (b) are illustrated in dot plots. Rich factor refers to the proportion of DMGs belonging to a specific term. Node size (gene number) refers to the number of DEGs within each term and node color indicates the level of significance (−log10 p value)

Additionally, we performed a pathway enrichment analysis using the Reactome database and DO analysis, and identified 117 and 55 significant pathways, respectively (Additional file 1: Figure S5). Enriched Reactome pathways primarily included those involved in insulin signaling such as “PI5P, PP2A and IER3 regulate PI3K/Akt signaling,” “PI3K cascade,” “IRS-mediated signaling,” and “IRS-related events triggered by IGF1R” (Additional file 1: Figure S5A, Additional file 2: Table S11). DMGs under Reactome were also associated with “VEGFA-VEGFR2 pathway” (Additional file 2: Table S11).

Overlap analysis between DMGs and DEGs

To understand the relationship between epigenetic regulation and transcriptomic changes in DPN and T2D, we examined the overlap between DMGs and DEGs. To that end, we screened 998 DEGs (Additional file 2: Table S3) and 1489 DMGs (Additional file 2: Table S8) for both altered gene expression and differential methylation, along with their corresponding directional changes. We found that 71 genes were shared between the DMGs and DEGs sets (Additional file 2: Table S13), and the majority of corresponding DMCs were located in intronic or intergenic regions. Among these 71 shared genes, 13 genes included multiple DMCs, while 58 genes included a single DMC per gene. These shared genes were divided into four categories based on the directional changes in DNA methylation and gene expression: “Hypo-Down” for the hypomethylated and downregulated genes (n = 20), “Hypo-Up” for the hypomethylated and upregulated genes (n = 33), “Hyper-Down” for the hypermethylated and downregulated genes (n = 14), and “Hyper-Up” for the hypermethylated and upregulated genes (n = 19) (Fig. 5a).

Fig. 5
figure5

Overlap between DMGs and DEGs. a Heatmap of methylation and expression changes of the overlapping DMGs and DEGs. The first column of each group corresponds to the methylation change (green: hypermethylated, blue: hypomethylated), while the second column represents the gene expression change (red: upregulated, blue: downregulated). DMCs belonging to the same genes are illustrated separately. b Dots include the significantly enriched KEGG pathways among the overlapping genes with opposite directions between DMGs and DEGs (hypo-up and hyper-down groups)

To further explore the interaction between DEGs and DMGs, we first constructed a gene interaction network using DEGs to identify gene sets that might highly interact with each other. DEGs were queried against the Search Tool for the Retrieval of Interacting Genes (STRING) database, and we generated an interaction network using the highest scoring confidence cutoff (score > 0.9). DMGs were then superimposed on the gene interaction network, and only the DEGs with a corresponding significant change in DNA methylation were listed (Fig. 6). This network included two large subnetworks and multiple small ones. The largest subnetwork included genes highly related to the immune response and locomotion and cell migration, and the second-largest subnetwork included genes involved in RNA binding and regulation of gene expression. Network analysis further showed a high degree of connectivity between phospholipase C gamma 2 (PLCG2) and G protein-coupled receptor 17 (GPR17), suggesting a role in DPN pathogenesis. Importantly, the network also highlighted highly connected overlapping genes between the differentially expressed and methylated gene sets, such as Thy-1 cell surface antigen (Thy1), also known as CD90. Thy1 has been implicated in nervous system development, neuronal injury, and immune response [37], and was found to be hypomethylated and downregulated in group 1 relative to group 2. Our results also show that mitogen-activated protein kinase 8 interacting protein 3 (MAPK8IP3), recently implicated in nerve degeneration and axonal neuropathy [38], is hypomethylated and upregulated in group 1 compared to group 2.

Fig. 6
figure6

Interaction network of DMGs and DEGs. An interaction network among the DEGs was generated using STRING, a database of known and predicted protein-protein interactions with the highest confidence score cutoff (score > 0.9). Node corresponds to a DEG, while edge indicates an interaction between two nodes. To visualize DMGs superimposed on the gene interaction network, only the DEGs with corresponding significant change in DNA methylation are labeled and highlighted with the color indicating the direction of methylation change: hypermethylation (red), hypomethylation (green), and both hyper- and hypomethylation (yellow). Major clusters were subjected to functional enrichment analysis in terms of GO and KEGG/Reactome pathways, with the most significantly enriched biological functions annotated in the network. Singletons not connected to another gene were excluded from this network

We lastly performed functional enrichment analysis on the genes showing opposite direction (Fig. 5b) and same direction (Additional file 1: Figure S6) of change in DNA methylation and gene expression. KEGG analysis showed that these genes are highly enriched in pathways including “ECM-receptor interaction,” “Pathways in cancer,” “PI3K−Akt signaling pathway,” and “calcium signaling pathway.” These numerous epigenetically regulated pathways are in agreement with our previous study, which demonstrated that many of the pathways previously identified as key mechanisms in human DPN are under epigenetic control [24].

Experimental validation using RT-qPCR

Based on biological relevance, we next selected two targets from the DMG and DEG overlap and functional enrichment analysis to validate by real-time quantitative PCR (RT-qPCR) (Fig. 7): Thy1 and prostaglandin-endoperoxide synthase 1 (PTGS1), with CpG sites located within intergenic regions. Thy1 was significantly hypomethylated by RRBS and downregulated by RNA-seq, and its gene expression followed a similar decreasing pattern in group 1 compared to group 2 (p = 0.0163; Fig. 7a). The directions of methylation and expression changes for PTGS1, a gene involved in neuropathic pain [39], were also discordant (hypomethylated/downregulated), and this finding was further confirmed by RT-qPCR, with a significant decrease in PTGS1 gene expression in group 1 versus group 2 (p = 0.0111; Fig. 7b).

Fig. 7
figure7

qPCR validation of genes with a similar direction of change in DNA methylation and gene expression. RNA was extracted from human sural nerve biopsies from group 1 (n  =  9–10) and group 2 (n  =  10), which was used as the relative control group (set to 100%), and quantified by RT-qPCR. Data were normalized to tyrosine 3-monooxygenase/tryptophan 5-monooxygenase activation protein zeta (YWHAZ) and expressed as % change calculated by the 2−ΔΔCT method. *p < 0.05 relative to group 2

Discussion

Although several genes and pathways are implicated in the pathophysiology of DPN, the molecular mechanisms underlying disease onset and progression remain largely unknown. DNA methylation, a major regulator of gene expression, has recently emerged as a key player in the development of diabetic complications, including DPN [24,25,26,27]. The aim of the present study was therefore to define the methylomic and transcriptomic signatures of human DPN and to understand how DNA methylation influences gene expression and thus contributes to nerve fiber damage. Using RNA-seq and RRBS, we determined differentially altered molecular pathways in both the methylome and transcriptome of sural nerves from well-characterized patients with T2D and DPN, and confirm findings gathered from experimental and clinical DPN [10, 24, 40, 41]. We also report the first integrated analysis of methylomic and transcriptomic datasets and identified differences in the methylation of genes encoding pathways involved in “immune response,” “ECM-receptor interaction,” and “PI3K−Akt signaling pathway.” These findings not only shed light on the role of epigenetic mechanisms in driving the expression of well-known regulators and novel targets of DPN but also help identify potential disease-modifying targets.

DPN is characterized by decreased nerve conduction velocities and the loss of both unmyelinated and myelinated nerve fibers [30, 40]. Accordingly, our initial analysis of the present cohort classified sural nerve biopsies based on their MFD, to reflect disease progression [30]. In this study, unbiased clustering of RNA-seq data separated those samples into distinctive clusters, which significantly differed in HbA1c levels. T2D patients with high HbA1c (8.7 ± 1.6%) showed distinct variations in their sural nerve transcriptome relative to patients with lower HbA1c (7.8 ± 1.4%). These findings suggest that transcriptomic changes in this cohort are associated with glycemic status, an established and independent risk factor for DPN development [30, 42].

When looking at transcriptomic changes, we found alterations in the expression of genes involved in immune response, calcium signaling, and axon guidance, which are highly relevant to nerve injury based on our current understanding of DPN pathogenesis [40, 43]. Moreover, the downregulation of the antioxidant SOD2 was of particular interest because this effect has been previously associated with increased oxidative damage and a stronger neuropathy phenotype in animal DPN models [44]. Consistent with these findings, the downregulated SOD2 suggests that antioxidant capacity is depleted in sural nerve biopsies of T2D patients with higher HbA1c and may participate in hyperglycemia-induced nerve injury. HbA1c-related nerve injury was also accompanied by changes in DNA methylation that mostly occurred within the gene body, a common pattern of diabetic complications [45]. DMGs were highly enriched for functions related to MAPK signaling, axon guidance, and VEGFA-VEGFR2 pathway. These observations, in line with our previous findings in a smaller human cohort and animal models of DPN, suggest a direct role for these pathways in DPN progression [10, 24].

Although RNA-seq and RRBS revealed transcriptomic and methylomic changes in sural nerve biopsies, this approach alone did not identify specific genes whose expression levels may be influenced by epigenetic factors. For the first time, we integrated methylation and gene expression datasets to determine whether the interaction between the methylome and transcriptome differed in T2D DPN patients with high (poor glycemic control) versus low HbA1c (good glycemic control). Specifically, our results determined that DNA methylation within the promoter or gene body was both concordantly and discordantly associated with gene expression. DNA methylation within promoter regions is associated with gene silencing and is generally considered a hallmark of cancer [21]. However, results from our group and others show that DNA methylation changes, in particular discordant changes, occur more frequently within gene bodies in diabetic complications, including DPN [24, 45]. Here, we confirmed and expanded our previous findings and demonstrated that the direction of change between DNA methylation and gene expression in human nerves can be both concordant and discordant.

Network analysis of overlapping DMG and DEG genes identified pathways involved in immune response, an important player in DPN development [46]. However, our results extend the literature to identify for the first time the role of DNA methylation in regulating immune-associated gene expression in human DPN. We found that PLCG2 was hypomethylated and overexpressed in sural nerves from patients with higher HbA1c. PLCG2 encodes a signaling protein essential for regulating immune cells, including macrophages [47], and has been recently implicated in neurodegenerative disorders [48, 49]. Peripheral nerve tissue consists of multiple cell types, including immune system macrophages and supportive Schwann cells, the myelinating cells of the peripheral nervous system, and axonal extensions. The current analysis did not account for the differential cell-type composition of sural nerve biopsies between groups 1 and 2. Therefore, it is unclear whether PLCG2 was expressed by the Schwann cells or infiltrating macrophages from the nerve biopsy tissue. However, evidence suggests that these two cell types interact with each other in the context of nerve injury and demyelination, and a similar mechanism may be occurring in human DPN in a PLCG2-associated manner [50, 51]. Future experiments addressing cell-specific changes would provide new insights into DPN pathogenesis.

Other immune-related genes of interest were GNAS and MAPK8IP3. GNAS, a complex imprinted locus with multiple gene products, including the G protein α-subunit Gsα [52], was hypermethylated and downregulated. GNAS is an interesting finding because it is critical for peripheral nerve myelination [53] and modulates lipid and glucose metabolism [54, 55], which are both dysregulated in experimental and clinical DPN [10, 14, 40]. We also detected hypomethylated and upregulated MAPK8IP3, a scaffold protein for c-Jun N-terminal kinase (JNK), in sural nerves of patients with higher HbA1c. In addition to its emerging role in nerve degeneration and axonal neuropathy [38], MAPK8IP3 may be an important player in the progression of human DPN through its close interaction with Toll-like receptor 4 (TLR4) and JNKs [56]. We recently demonstrated a role for TLR4 in immune-mediated inflammation in murine DPN [57]. Additionally, JNK is a key signaling pathway promoting inflammation and insulin resistance in diabetic complications, including DPN [40, 58]. Thus, our results support a role for epigenetic mechanisms in regulating immune-associated genes such as PLCG2, GNAS, and MAPK8IP3 and suggest that they may be potential disease-modifying therapeutic targets in human DPN.

Another mechanism thought to contribute to DPN is ECM remodeling [59], and we have previously reported transcriptomic changes in the composition and function of the ECM pathways in DPN [9, 14]. In the current study, the overlapping genes between DMGs and DEGs were enriched with biological functions related to locomotion/cell migration and ECM-receptor interaction, consistent with the previous findings [9, 14]. In particular, we showed that the hypomethylation of tissue inhibitor of metalloproteinase 2 (TIMP2) was associated with reduced gene expression. TIMP2 regulates the activity of matrix metalloproteinase 2 (MMP2), and the MMP2/TIMP2 ratio is essential for maintaining ECM integrity [60]. An imbalance between MMP2 and TIMP2 contributes to ECM accumulation and fibrosis in diabetic nephropathy [61], an effect that may be mediated at least in part by TGF-β [62], which was also dysregulated in the present study. Additionally, the MMP2/TIMP2 axis modulates sciatic nerve ECM during nerve repair [63] and an impaired MMP2/TIMP2 axis has been implicated in the pathogenesis of DPN [64]. Given its prominent role in the regulation of nerve ECM, it is not surprising that our gene interaction network revealed that TIMP2 was highly connected to members of the collagen family, such as collagen type VI, α3 (COL6A3), which has also been associated with fibrosis and inflammation in diabetic nephropathy [65]. While we speculate that alteration of TIMP2 and COL6A3 in human DPN results in ECM remodeling, further investigation will determine how these changes impact nerve function in diabetes.

Impaired insulin signaling in DPN may induce myelination deficits in Schwann cells and insulin resistance in sensory neurons [66, 67]. Our results show that insulin signaling pathway genes such as IRS2 are differentially methylated and that IGF-I is downregulated, indicating impaired insulin signaling in the sural nerves of patients with higher HbA1c. Insulin and IGF-I both signal through the PI3K-Akt pathway, which in turn exerts multiple cellular actions through downstream effectors, including the mammalian target of rapamycin complex 1 (mTORC1) [68]. PI3K-Akt-mTORC1 signaling is heavily implicated in Schwann cell lipid synthesis, a critical mechanism for myelination [66], whose disruption is associated with impaired nerve function [69]. Consistent with these findings, our KEGG-based analysis revealed a dysregulated PI3K-Akt pathway at both the DNA methylation and gene expression levels, and suggests that DNA methylation may be a new mechanism for regulating PI3K-Akt in peripheral nerves. Although the mTORC1 pathway is not differentially expressed in this cohort, preliminary evidence from our group supports the involvement of mTORC1 in T2D-mediated nerve injury (data not shown), and future studies will elucidate the interaction between PI3K-Akt and mTORC1 in the context of DPN.

Thy1 and PTGS1 were the genes selected for RTq-PCR validation because of their relevance to known nerve injury mechanisms in DPN [37, 39]. Thy1 and PTGS1 were hypomethylated within the intergenic region and downregulated in gene expression, which is expected when methylation occurs outside the promoter region [22]. Chen et al. have shown that Thy1 levels are reduced following nerve crush and return to near control levels after nerve repair [70]. While the type of nerve insult is different in our study (metabolic versus acute), our results also show that Thy1 is downregulated following injury, supporting the idea that this gene is a negative regulator of neurite outgrowth. PTGS1 encodes cyclooxygenase 1 (COX-1), a modulator of inflammation, which influences neuropathic pain in dorsal root ganglia and spinal neurons [39, 71]. Inhibiting COX-1 in streptozotocin-induced diabetic rats attenuates hyperalgesia [72]. The hypomethylation-dependent reduction in PTGS1 gene expression in the sural nerves of patients with higher HbA1c suggests it is involved in DPN pathogenesis, supports our findings of a role for inflammation in DPN, and warrants further investigation.

Conclusions

We profiled the sural nerve methylome and transcriptome in human DPN to identify changes in molecular mechanisms underlying disease pathogenesis. We show that HbA1c levels are associated with transcriptomic changes and pathways which modulate important molecular functions, such as the immune response and axon guidance. Our functional and network analyses of integrated epigenetic and transcriptomic data revealed regulation of some of these pathways by DNA methylation, a potentially reversible mechanism linking genetics and lifestyle factors. We found that immune response, ECM-receptor interaction, and PI3K-Akt signaling pathways are under epigenetic control and may play a crucial role in DPN development. Although the exact mechanisms are yet to be elucidated, our results suggest that optimal glycemic control is one of the important factors in maintaining epigenetic homeostasis and nerve function.

Methods

Sample collection and preparation

Human sural nerve biopsies were collected during a previous 52-week double-blind placebo-controlled clinical trial of acetyl-L-carnitine for treating DPN [28, 29]. The trial design and selection criteria were described previously [29]. Briefly, the trial included adult male and female subjects diagnosed with T1D or T2D for at least 1 year and with an HbA1c > 5.9%, with clinically evident DPN, defined by at least 2 neurological findings, including clinical symptoms or abnormal electrophysiological tests (nerve conduction velocity or vibration perception threshold) [73, 74]. Patients with non-diabetic causes of neuropathy, complicating diseases (such as HIV or significant cardiac or hepatic disorders), and alcohol or drug abuse were excluded.

A sural nerve biopsy was collected from the distal calf at the time patients enrolled and another biopsy was collected from the opposite leg after 52 weeks of treatment. Sural nerve myelinated fiber density (MFD) is a reliable morphological criterion for DPN diagnosis [75]. Based on percent myelinated fiber density change (%delta-MFD) over 1 year, these patients were divided into three groups (denoted as regenerator, intermediator, and degenerator) in the previous analysis [30]. For the current study, a total of 78 patients with T2D were selected for methylome and transcriptome profiling; patient characteristics are given in Additional file 2: Table S1. Different from previous analyse s[30], these 78 patients were segregated into groups 1 and 2 using an unbiased clustering analysis (see below methods). The University of Michigan Institutional Review Board for Human Subject Research approved the use of human sural nerve samples.

Genome-wide gene expression and methylation profiling

DNA and RNA were extracted from human sural nerves using the Qiagen AllPrep DNA/RNA Mini Kit (Qiagen, Valencia, CA, USA) according to the manufacturer’s protocol. DNA and RNA quantity was assessed using a Qubit fluorometer, and the quality was evaluated using a 2100 Bioanalyzer (Agilent Technologies, Santa Clara, CA).

Samples with RNA integrity numbers ≥ 8 were prepared for RNA-sequencing (RNA-seq) using the Illumina TruSeq mRNA Sample Prep v2 kit (Illumina, San Diego, CA, USA). Approximately 90 million 50-base pair (bp) paired-end reads per sample were obtained on an Illumina HiSeq 2500 system (Illumina, San Diego, CA). RNA-seq was performed by the University of Michigan DNA Sequencing Core (http://seqcore.brcf.ed.umich.edu).

RRBS was performed by the University of Michigan Epigenomics Core as previously described [76]. Briefly, DNA was digested with the MspI restriction enzyme and purified using phenol:chloroform extraction and ethanol precipitation, blunt-end digested, phosphorylated, and ligated into an adapter duplex with methylated cytosines. The ligated fragments were cleaned and size selected by an agarose gel. Bisulfite conversion was performed for selected fragments followed by PCR amplification and cleanup using AMPure XP beads. The Qubit assay and Agilent’s High Sensitivity D1000 assay were used to quantify the libraries, which were then sequenced on the Illumina HiSeq 2500 platform to obtain approximately 90 million 50 bp single-end reads per sample.

Sequencing data analysis

Quality filtering and read mapping

For the RNA-seq data, low-quality (Q < 30) reads and sequencing adapters were removed with Trimmomatic [77] and the quality of raw reads was assessed with FastQC (version 0.11.5, Babraham Bioinformatics, UK). Clean reads were then mapped to the human reference genome hg19 Refseq using Hisat2 [78]. FeatureCounts [79] was used to count reads mapped to individual genes and only uniquely mapped reads were used in the counting step. Genes with zero expression across all samples were omitted from the correlation and differential expression analysis. Fragments per kilobase of exon per million mapped reads (FPKM) was calculated to represent gene expression levels.

Similarly for the RRBS data, quality control of the RRBS data was performed using FastQC (version 0.11.5, Babraham Bioinformatics, UK), and low-quality reads were removed with Trim Galore (version 0.5.0, Babraham Bioinformatics, UK). Then, trimmed reads were aligned and mapped using Bismark [80] to the human hg19 reference genome. The percentage methylation level was calculated by #C/(#C + #T), where #C is the number of methylated reads and #T is the unmethylated reads. Then, only the CpG sites with a read coverage > 10, a quality score > 20, and appeared at least in 10 samples among each group, using the parameter “min.per.group = 10,” were kept for downstream analysis.

Unbiased sample clustering based on RNA-seq data

The RNA-seq data was normalized by DESeq2 R package using the default parameters [81]. Unsupervised hierarchical clustering and principal component analysis on the normalized expression values were used to examine the overall similarity among the samples, which identified three sample groups with high similarity. We also determined whether clinical factors (Additional file 2: Table S1) were significantly associated with the three groups using multifactorial logistic regression analyses.

Differential expression and methylation analyses

Differentially expressed genes (DEGs) were identified using DESeq2 [81] and genes with a Benjamini-Hochberg adjusted p value < 0.01 were deemed as DEGs. Differentially methylated CpGs (DMCs) were defined as a CpG site with a methylation difference of > 15% and a false discovery rate adjusted p value (q value) < 0.01 between group 1 and group 2. DMCs were then annotated based on genes and CpG island (CGi) features. Each DMC was mapped to a gene, having the shortest distance from its transcript starting site to the DMC. These mapped genes were defined as differentially methylated genes (DMGs).

Functional enrichment analysis

To identify and compare the overrepresented biological functions, enrichment analysis was performed using a hypergeometric test with our in-house R analysis package richR (http://github.com/hurlab/richR). Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways, Gene Ontology (GO) terms, Reactome pathway (https://reactome.org), and Disease Ontology (DO: http://disease-ontology.org) terms were used in the enrichment analysis, and the terms with Benjamini-Hochberg corrected p values < 0.05 were deemed as significantly overrepresented biological functions in each DEG set. To reduce redundancy in the GO enrichment results, additional term clustering was performed using the clustering parameter kappa > 0.5, which included GO clusters with a minimum of 5 GO terms, as implemented in the Database for Annotation, Visualization and Integrated Discovery (DAVID) [33, 82]. The same functional enrichment analysis was performed for the DMG sets with a nominal p value < 0.05 as the cutoff value. Dot plots with the top 20 most significant terms were generated.

RT-qPCR validation

Two shared genes between DMGs and DEGs (Thy1 and PTGS1) were selected for validation by quantitative real-time polymerase chain reaction (RT-qPCR). RNA was isolated from sural nerve biopsies using RNeasy fibrous tissue mini kit (Qiagen, Valencia, CA, USA) from group 1 (n  =  9–10) and group 2 (n  =  10), which was used as the relative control group (set to 100%). Reverse transcription was performed using iScript Supermix (Bio-Rad Laboratories, Hercules, California). qPCR reactions were carried out using sequence-specific TaqMan™ probes (ThermoFisher/Applied Biosystems) for Thy1 (Hs06633377_s1) and PTGS1 (Hs00377726_m1) in an Applied Biosystems StepOneTM RT-PCR system. Using the 2−ΔΔCT method, tyrosine 3-monooxygenase/tryptophan 5-monooxygenase activation protein (YWHAZ) was used as the endogenous control and group 2 as the relative control. Statistically significant differences were calculated using a Student’s t test with the GraphPad Prism 7 software (GraphPad Software Inc.).

Availability of data and materials

The datasets generated and/or analyzed during the current study are available in the NCBI Gene Expression Omnibus (Accession ID: GSE148061; token for reviewer’s access: “cfytwymepnwjdqf” without quotes).

Change history

  • 27 August 2020

    An amendment to this paper has been published and can be accessed via the original article.

Abbreviations

DAVID:

Database for Annotation, Visualization and Integrated Discovery

DEG:

Differentially expressed gene

DMC:

Differentially methylated CpG site

DMG:

Differentially methylated gene

DO:

Disease Ontology

DPN:

Diabetic peripheral neuropathy

GO:

Gene Ontology

KEGG:

Kyoto Encyclopedia of Genes and Genomes

MFD:

Myelinated fiber density

RNA-seq:

RNA-sequencing

RRBS:

Reduced representation bisulfite sequencing

STRING:

Search Tool for the Retrieval of Interacting Genes

TSS:

Transcription starting site

References

  1. 1.

    Saeedi P, Petersohn I, Salpea P, Malanda B, Karuranga S, Unwin N, Colagiuri S, Guariguata L, Motala AA, Ogurtsova K, et al: Global and regional diabetes prevalence estimates for 2019 and projections for 2030 and 2045: results from the International Diabetes Federation Diabetes Atlas, 9(th) edition. Diabetes Res Clin Pract 2019, 157:107843.

  2. 2.

    Zheng Y, Ley SH, Hu FB. Global aetiology and epidemiology of type 2 diabetes mellitus and its complications. Nat Rev Endocrinol. 2018;14:88–98.

    PubMed  PubMed Central  Google Scholar 

  3. 3.

    Chatterjee S, Khunti K, Davies MJ. Type 2 diabetes. Lancet. 2017;389:2239–51.

    CAS  PubMed  PubMed Central  Google Scholar 

  4. 4.

    Callaghan BC, Price RS, Feldman EL. Distal symmetric polyneuropathy: a review. JAMA. 2015;314:2172–81.

    CAS  PubMed  PubMed Central  Google Scholar 

  5. 5.

    Carinci F, Massi Benedetti M, Klazinga NS, Uccioli L. Lower extremity amputation rates in people with diabetes as an indicator of health systems performance. A critical appraisal of the data collection 2000-2011 by the Organization for Economic Cooperation and Development (OECD). Acta Diabetol. 2016;53:825–32.

    CAS  PubMed  PubMed Central  Google Scholar 

  6. 6.

    Edwards JL, Vincent AM, Cheng HT, Feldman EL. Diabetic neuropathy: mechanisms to management. Pharmacol Ther. 2008;120:1–34.

    CAS  PubMed  PubMed Central  Google Scholar 

  7. 7.

    Callaghan BC, Little AA, Feldman EL, Hughes RA. Enhanced glucose control for preventing and treating diabetic neuropathy. Cochrane Database Syst Rev. 2012:CD007543.

  8. 8.

    Papanas N, Ziegler D. Risk factors and comorbidities in diabetic neuropathy: an update 2015. Rev Diabet Stud. 2015;12:48–62.

    PubMed  PubMed Central  Google Scholar 

  9. 9.

    Hinder LM, Murdock BJ, Park M, Bender DE, O'Brien PD, Rumora AE, Hur J, Feldman EL. Transcriptional networks of progressive diabetic peripheral neuropathy in the db/db mouse model of type 2 diabetes: an inflammatory story. Exp Neurol. 2018;305:33–43.

    CAS  PubMed  PubMed Central  Google Scholar 

  10. 10.

    Pande M, Hur J, Hong Y, Backus C, Hayes JM, Oh SS, Kretzler M, Feldman EL. Transcriptional profiling of diabetic neuropathy in the BKS db/db mouse: a model of type 2 diabetes. Diabetes. 2011;60:1981–9.

    CAS  PubMed  PubMed Central  Google Scholar 

  11. 11.

    Hinder LM, Park M, Rumora AE, Hur J, Eichinger F, Pennathur S, Kretzler M, Brosius FC 3rd, Feldman EL. Comparative RNA-Seq transcriptome analyses reveal distinct metabolic pathways in diabetic nerve and kidney disease. J Cell Mol Med. 2017;21:2140–52.

    CAS  PubMed  PubMed Central  Google Scholar 

  12. 12.

    Ma J, Pan P, Anyika M, Blagg BS, Dobrowsky RT. Modulating molecular chaperones improves mitochondrial bioenergetics and decreases the inflammatory transcriptome in diabetic sensory neurons. ACS Chem Neurosci. 2015;6:1637–48.

    CAS  PubMed  PubMed Central  Google Scholar 

  13. 13.

    Hur J, Dauch JR, Hinder LM, Hayes JM, Backus C, Pennathur S, Kretzler M, Brosius FC 3rd, Feldman EL. The metabolic syndrome and microvascular complications in a murine model of type 2 diabetes. Diabetes. 2015;64:3294–304.

    CAS  PubMed  PubMed Central  Google Scholar 

  14. 14.

    McGregor BA, Eid S, Rumora AE, Murdock B, Guo K, de Anda-Jauregui G, Porter JE, Feldman EL, Hur J. Conserved transcriptional signatures in human and murine diabetic peripheral neuropathy. Sci Rep. 2018;8:17678.

    PubMed  PubMed Central  Google Scholar 

  15. 15.

    Sas KM, Lin J, Rajendiran TM, Soni T, Nair V, Hinder LM, Jagadish HV, Gardner TW, Abcouwer SF, Brosius FC 3rd, et al. Shared and distinct lipid-lipid interactions in plasma and affected tissues in a diabetic mouse model. J Lipid Res. 2018;59:173–83.

    CAS  PubMed  PubMed Central  Google Scholar 

  16. 16.

    Sas KM, Kayampilly P, Byun J, Nair V, Hinder LM, Hur J, Zhang H, Lin C, Qi NR, Michailidis G, et al. Tissue-specific metabolic reprogramming drives nutrient flux in diabetic complications. JCI Insight. 2016;1:e86976.

    PubMed  PubMed Central  Google Scholar 

  17. 17.

    Mantione KJ, Kream RM, Kuzelova H, Ptacek R, Raboch J, Samuel JM, Stefano GB. Comparing bioinformatic gene expression profiling methods: microarray and RNA-Seq. Med Sci Monit Basic Res. 2014;20:138–42.

    PubMed  PubMed Central  Google Scholar 

  18. 18.

    Malone JH, Oliver B. Microarrays, deep sequencing and the true measure of the transcriptome. BMC Biol. 2011;9:34.

    CAS  PubMed  PubMed Central  Google Scholar 

  19. 19.

    Quackenbush J. Microarray data normalization and transformation. Nat Genet. 2002;32(Suppl):496–501.

    CAS  PubMed  PubMed Central  Google Scholar 

  20. 20.

    Everaert C, Luypaert M, Maag JLV, Cheng QX, Dinger ME, Hellemans J, Mestdagh P. Benchmarking of RNA-sequencing analysis workflows using whole-transcriptome RT-qPCR expression data. Sci Rep. 2017;7:1559.

    PubMed  PubMed Central  Google Scholar 

  21. 21.

    Maunakea AK, Chepelev I, Zhao K. Epigenome mapping in normal and disease States. Circ Res. 2010;107:327–39.

    CAS  PubMed  PubMed Central  Google Scholar 

  22. 22.

    Greenberg MVC, Bourc'his D. The diverse roles of DNA methylation in mammalian development and disease. Nat Rev Mol Cell Biol. 2019;20:590–607.

    CAS  PubMed  PubMed Central  Google Scholar 

  23. 23.

    Jones PA. Functions of DNA methylation: islands, start sites, gene bodies and beyond. Nat Rev Genet. 2012;13:484–92.

    CAS  PubMed  PubMed Central  Google Scholar 

  24. 24.

    Guo K, Elzinga S, Eid S, Figueroa-Romero C, Hinder LM, Pacut C, Feldman EL, Hur J. Genome-wide DNA methylation profiling of human diabetic peripheral neuropathy in subjects with type 2 diabetes mellitus. Epigenetics. 2019;14:766–79.

    PubMed  PubMed Central  Google Scholar 

  25. 25.

    Davegardh C, Garcia-Calzon S, Bacos K, Ling C. DNA methylation in the pathogenesis of type 2 diabetes in humans. Mol Metab. 2018;14:12–25.

    PubMed  PubMed Central  Google Scholar 

  26. 26.

    VanderJagt TA, Neugebauer MH, Morgan M, Bowden DW, Shah VO. Epigenetic profiles of pre-diabetes transitioning to type 2 diabetes and nephropathy. World J Diabetes. 2015;6:1113–21.

    PubMed  PubMed Central  Google Scholar 

  27. 27.

    Maghbooli Z, Hossein-nezhad A, Larijani B, Amini M, Keshtkar A. Global DNA methylation as a possible biomarker for diabetic retinopathy. Diabetes Metab Res Rev. 2015;31:183–9.

    CAS  PubMed  PubMed Central  Google Scholar 

  28. 28.

    Wiggin TD, Sullivan KA, Pop-Busui R, Amato A, Sima AA, Feldman EL. Elevated triglycerides correlate with progression of diabetic neuropathy. Diabetes. 2009;58:1634–40.

    CAS  PubMed  PubMed Central  Google Scholar 

  29. 29.

    Sima AA, Calvani M, Mehra M, Amato A. Acetyl LCSG: acetyl-L-carnitine improves pain, nerve regeneration, and vibratory perception in patients with chronic diabetic neuropathy: an analysis of two randomized placebo-controlled trials. Diabetes Care. 2005;28:89–94.

    CAS  PubMed  PubMed Central  Google Scholar 

  30. 30.

    Hur J, Sullivan KA, Callaghan BC, Pop-Busui R, Feldman EL. Identification of factors associated with sural nerve regeneration and degeneration in diabetic neuropathy. Diabetes Care. 2013;36:4043–9.

    CAS  PubMed  PubMed Central  Google Scholar 

  31. 31.

    Prinz F, Kapeller A, Pichler M, Klec C. The implications of the long non-coding RNA NEAT1 in non-cancerous diseases. Int J Mol Sci. 2019;20.

  32. 32.

    Lin Z, Li X, Zhan X, Sun L, Gao J, Cao Y, Qiu H. Construction of competitive endogenous RNA network reveals regulatory role of long non-coding RNAs in type 2 diabetes mellitus. J Cell Mol Med. 2017;21:3204–13.

    CAS  PubMed  PubMed Central  Google Scholar 

  33. 33.

    Huang da W, Sherman BT, Lempicki RA: Systematic and integrative analysis of large gene lists using DAVID bioinformatics resources. Nat Protoc 2009, 4:44-57.

  34. 34.

    Anjaneyulu M, Berent-Spillson A, Inoue T, Choi J, Cherian K, Russell JW. Transforming growth factor-beta induces cellular injury in experimental diabetic neuropathy. Exp Neurol. 2008;211:469–79.

    CAS  PubMed  PubMed Central  Google Scholar 

  35. 35.

    Flynn JM, Melov S. SOD2 in mitochondrial dysfunction and neurodegeneration. Free Radic Biol Med. 2013;62:4–12.

    CAS  PubMed  PubMed Central  Google Scholar 

  36. 36.

    Sullivan KA, Kim B, Feldman EL. Insulin-like growth factors in the peripheral nervous system. Endocrinology. 2008;149:5963–71.

    CAS  PubMed  PubMed Central  Google Scholar 

  37. 37.

    Bradley JE, Ramirez G, Hagood JS. Roles and regulation of Thy-1, a context-dependent modulator of cell phenotype. Biofactors. 2009;35:258–65.

    CAS  PubMed  PubMed Central  Google Scholar 

  38. 38.

    Iwasawa SYK, Kikuchi A, Kobayashi Y, Haginoya K, Matsumoto H, Kurosawa K, Ochiai M, Sakai Y, Fujita A, et al. Recurrent de novo MAPK8IP3 variants cause neurological phenotypes. Ann Neurol. 2019;85:927–33.

    CAS  PubMed  PubMed Central  Google Scholar 

  39. 39.

    Kanda H, Kobayashi K, Yamanaka H, Noguchi K. COX-1-dependent prostaglandin D2 in microglia contributes to neuropathic pain via DP2 receptor in spinal neurons. Glia. 2013;61:943–56.

    PubMed  PubMed Central  Google Scholar 

  40. 40.

    Feldman EL, Nave KA, Jensen TS, Bennett DLH. New horizons in diabetic neuropathy: mechanisms, bioenergetics, and pain. Neuron. 2017;93:1296–313.

    CAS  PubMed  PubMed Central  Google Scholar 

  41. 41.

    Hur JSK, Pande M, Hong Y, Sima AA, Jagadish HV, Kretzler M, Feldman EL. The identification of gene expression profiles associated with progression of human diabetic neuropathy. Brain. 2011;134:3222–35.

    PubMed  PubMed Central  Google Scholar 

  42. 42.

    Callaghan BCGL, Li Y, Zhou X, Reynolds E, Banerjee M, Pop-Busui R, Feldman EL, Ji L. Diabetes and obesity are the main metabolic drivers of peripheral neuropathy. Ann Clin Transl Neurol. 2018;5:397–405.

    PubMed  PubMed Central  Google Scholar 

  43. 43.

    Fernyhough P, Calcutt NA. Abnormal calcium homeostasis in peripheral neuropathies. Cell Calcium. 2010;47:130–9.

    CAS  PubMed  PubMed Central  Google Scholar 

  44. 44.

    Vincent AM, Russell JW, Sullivan KA, Backus C, Hayes JM, McLean LL, Feldman EL. SOD2 protects neurons from injury in cell culture and animal models of diabetic neuropathy. Exp Neurol. 2007;208:216–27.

    CAS  PubMed  PubMed Central  Google Scholar 

  45. 45.

    Chen Z, Miao F, Paterson AD, Lachin JM, Zhang L, Schones DE, Wu X, Wang J, Tompkins JD, Genuth S, et al. Epigenomic profiling reveals an association between persistence of DNA methylation and metabolic memory in the DCCT/EDIC type 1 diabetes cohort. Proc Natl Acad Sci U S A. 2016;113:E3002–11.

    CAS  PubMed  PubMed Central  Google Scholar 

  46. 46.

    Pop-Busui R, Ang L, Holmes C, Gallagher K, Feldman EL. Inflammation as a therapeutic target for diabetic neuropathies. Curr Diab Rep. 2016;16:29.

    PubMed  PubMed Central  Google Scholar 

  47. 47.

    Di Raimo T, Leopizzi M, Mangino G, Rocca CD, Businaro R, Longo L, Lo Vasco VR. Different expression and subcellular localization of phosphoinositide-specific phospholipase C enzymes in differently polarized macrophages. J Cell Commun Signal. 2016;10:283–93.

    PubMed  PubMed Central  Google Scholar 

  48. 48.

    van der Lee SJ, Conway OJ, Jansen I, Carrasquillo MM, Kleineidam L, van den Akker E, Hernandez I, van Eijk KR, Stringa N, Chen JA, et al. A nonsynonymous mutation in PLCG2 reduces the risk of Alzheimer’s disease, dementia with Lewy bodies and frontotemporal dementia, and increases the likelihood of longevity. Acta Neuropathol. 2019;138:237–50.

    PubMed  PubMed Central  Google Scholar 

  49. 49.

    Sims R, van der Lee SJ, Naj AC, Bellenguez C, Badarinarayan N, Jakobsdottir J, Kunkle BW, Boland A, Raybould R, Bis JC, et al. Rare coding variants in PLCG2, ABI3, and TREM2 implicate microglial-mediated innate immunity in Alzheimer’s disease. Nat Genet. 2017;49:1373–84.

    CAS  PubMed  PubMed Central  Google Scholar 

  50. 50.

    Stratton JA, Holmes A, Rosin NL, Sinha S, Vohra M, Burma NE, Trang T, Midha R, Biernaskie J. Macrophages regulate Schwann cell maturation after nerve injury. Cell Rep. 2018;24:2561–72 e2566.

    CAS  PubMed  PubMed Central  Google Scholar 

  51. 51.

    Martini R, Fischer S, Lopez-Vales R, David S. Interactions between Schwann cells and macrophages in injury and inherited demyelinating disease. Glia. 2008;56:1566–77.

    PubMed  PubMed Central  Google Scholar 

  52. 52.

    Weinstein LS, Xie T, Qasem A, Wang J, Chen M. The role of GNAS and other imprinted genes in the development of obesity. Int J Obes (Lond). 2010;34:6–17.

    CAS  Google Scholar 

  53. 53.

    Deng Y, Wu LMN, Bai S, Zhao C, Wang H, Wang J, Xu L, Sakabe M, Zhou W, Xin M, Lu QR. A reciprocal regulatory loop between TAZ/YAP and G-protein Galphas regulates Schwann cell proliferation and myelination. Nat Commun. 2017;8:15161.

    PubMed  PubMed Central  Google Scholar 

  54. 54.

    Chen M, Haluzik M, Wolf NJ, Lorenzo J, Dietz KR, Reitman ML, Weinstein LS. Increased insulin sensitivity in paternal Gnas knockout mice is associated with increased lipid clearance. Endocrinology. 2004;145:4094–102.

    CAS  PubMed  PubMed Central  Google Scholar 

  55. 55.

    Chen M, Gavrilova O, Liu J, Xie T, Deng C, Nguyen AT, Nackers LM, Lorenzo J, Shen L, Weinstein LS. Alternative Gnas gene products have opposite effects on glucose and lipid metabolism. Proc Natl Acad Sci U S A. 2005;102:7386–91.

    CAS  PubMed  PubMed Central  Google Scholar 

  56. 56.

    Matsuguchi T, Masuda A, Sugimoto K, Nagai Y, Yoshikai Y. JNK-interacting protein 3 associates with Toll-like receptor 4 and is involved in LPS-mediated JNK activation. EMBO J. 2003;22:4455–64.

    CAS  PubMed  PubMed Central  Google Scholar 

  57. 57.

    Elzinga S, Murdock BJ, Guo K, Hayes JM, Tabbey MA, Hur J, Feldman EL. Toll-like receptors and inflammation in metabolic neuropathy; a role in early versus late disease? Exp Neurol. 2019;320:112967.

    CAS  PubMed  PubMed Central  Google Scholar 

  58. 58.

    Solinas G, Becattini B. JNK at the crossroad of obesity, insulin resistance, and cell stress response. Mol Metab. 2017;6:174–84.

    CAS  PubMed  PubMed Central  Google Scholar 

  59. 59.

    Hill R. Extracellular matrix remodelling in human diabetic neuropathy. J Anat. 2009;214:219–25.

    PubMed  PubMed Central  Google Scholar 

  60. 60.

    Bernardo MM, Fridman R. TIMP-2 (tissue inhibitor of metalloproteinase-2) regulates MMP-2 (matrix metalloproteinase-2) activity in the extracellular environment after pro-MMP-2 activation by MT1 (membrane type 1)-MMP. Biochem J. 2003;374:739–45.

    CAS  PubMed  PubMed Central  Google Scholar 

  61. 61.

    Garcia-Fernandez N, Jacobs-Cacha C, Mora-Gutierrez JM, Vergara A, Orbe J, Soler MJ. Matrix metalloproteinases in diabetic kidney disease. J Clin Med. 2020;9.

  62. 62.

    Han SY, Jee YH, Han KH, Kang YS, Kim HK, Han JY, Kim YS, Cha DR. An imbalance between matrix metalloproteinase-2 and tissue inhibitor of matrix metalloproteinase-2 contributes to the development of early diabetic nephropathy. Nephrol Dial Transplant. 2006;21:2406–16.

    CAS  PubMed  PubMed Central  Google Scholar 

  63. 63.

    Gantus MA, Nasciutti LE, Cruz CM, Persechini PM, Martinez AM. Modulation of extracellular matrix components by metalloproteinases and their tissue inhibitors during degeneration and regeneration of rat sural nerve. Brain Res. 2006;1122:36–46.

    CAS  PubMed  PubMed Central  Google Scholar 

  64. 64.

    Ali S, Driscoll HE, Newton VL, Gardiner NJ. Matrix metalloproteinase-2 is downregulated in sciatic nerve by streptozotocin induced diabetes and/or treatment with minocycline: implications for nerve regeneration. Exp Neurol. 2014;261:654–65.

    CAS  PubMed  PubMed Central  Google Scholar 

  65. 65.

    Woroniecka KI, Park AS, Mohtat D, Thomas DB, Pullman JM, Susztak K. Transcriptome analysis of human diabetic kidney disease. Diabetes. 2011;60:2354–69.

    CAS  PubMed  PubMed Central  Google Scholar 

  66. 66.

    Hackett AR, Strickland A, Milbrandt J. Disrupting insulin signaling in Schwann cells impairs myelination and induces a sensory neuropathy. Glia. 2019.

  67. 67.

    Grote CW, Morris JK, Ryals JM, Geiger PC, Wright DE. Insulin receptor substrate 2 expression and involvement in neuronal insulin resistance in diabetic neuropathy. Exp Diabetes Res. 2011;2011:212571.

    CAS  PubMed  PubMed Central  Google Scholar 

  68. 68.

    Kim B, Feldman EL. Insulin resistance in the nervous system. Trends Endocrinol Metab. 2012;23:133–41.

    CAS  PubMed  PubMed Central  Google Scholar 

  69. 69.

    Sherman DL, Krols M, Wu LM, Grove M, Nave KA, Gangloff YG, Brophy PJ. Arrest of myelination and reduced axon growth when Schwann cells lack mTOR. J Neurosci. 2012;32:1817–25.

    CAS  PubMed  PubMed Central  Google Scholar 

  70. 70.

    Chen CH, Wang SM, Yang SH, Jeng CJ. Role of Thy-1 in in vivo and in vitro neural development and regeneration of dorsal root ganglionic neurons. J Cell Biochem. 2005;94:684–94.

    CAS  PubMed  PubMed Central  Google Scholar 

  71. 71.

    Dou W, Jiao Y, Goorha S, Raghow R, Ballou LR. Nociception and the differential expression of cyclooxygenase-1 (COX-1), the COX-1 variant retaining intron-1 (COX-1v), and COX-2 in mouse dorsal root ganglia (DRG). Prostaglandins Other Lipid Mediat. 2004;74:29–43.

    CAS  PubMed  PubMed Central  Google Scholar 

  72. 72.

    Bujalska M, Tatarkiewicz J, de Corde A, Gumulka SW. Effect of cyclooxygenase and nitric oxide synthase inhibitors on streptozotocin-induced hyperalgesia in rats. Pharmacology. 2008;81:151–7.

    CAS  PubMed  PubMed Central  Google Scholar 

  73. 73.

    Greene DA, Arezzo JC, Brown MB. Effect of aldose reductase inhibition on nerve conduction and morphometry in diabetic neuropathy. Zenarestat Study Group. Neurology. 1999;53:580–91.

  74. 74.

    O'Brien PC. Procedures for comparing samples with multiple endpoints. Biometrics. 1984;40:1079–87.

    CAS  PubMed  PubMed Central  Google Scholar 

  75. 75.

    Behse F, Buchthal F, Carlsen F. Nerve biopsy and conduction studies in diabetic neuropathy. J Neurol Neurosurg Psychiatry. 1977;40:1072–82.

    CAS  PubMed  PubMed Central  Google Scholar 

  76. 76.

    Day SE, Coletta RL, Kim JY, Campbell LE, Benjamin TR, Roust LR, De Filippis EA, Dinu V, Shaibi GQ, Mandarino LJ, Coletta DK. Next-generation sequencing methylation profiling of subjects with obesity identifies novel gene changes. Clin Epigenetics. 2016;8:77.

    PubMed  PubMed Central  Google Scholar 

  77. 77.

    Bolger AM, Lohse M, Usadel B. Trimmomatic: a flexible trimmer for Illumina sequence data. Bioinformatics. 2014;30:2114–20.

    CAS  PubMed  PubMed Central  Google Scholar 

  78. 78.

    Kim D, Langmead B, Salzberg SL. HISAT: a fast spliced aligner with low memory requirements. Nat Methods. 2015;12:357–60.

    CAS  PubMed  PubMed Central  Google Scholar 

  79. 79.

    Liao Y, Smyth GK. Shi W: featureCounts: an efficient general purpose program for assigning sequence reads to genomic features. Bioinformatics. 2014;30:923–30.

    CAS  PubMed  PubMed Central  Google Scholar 

  80. 80.

    Krueger F, Andrews SR. Bismark: a flexible aligner and methylation caller for Bisulfite-Seq applications. Bioinformatics. 2011;27:1571–2.

    CAS  PubMed  PubMed Central  Google Scholar 

  81. 81.

    Love MI, Huber W, Anders S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol. 2014;15:550.

    PubMed  PubMed Central  Google Scholar 

  82. 82.

    Huang da W, Sherman BT, Lempicki RA: Bioinformatics enrichment tools: paths toward the comprehensive functional analysis of large gene lists. Nucleic Acids Res 2009, 37:1-13.

Download references

Acknowledgements

The authors would like to thank Drs. Stacey Jacoby and Masha Savieleff for expert assistance in data interpretation and manuscript editing.

Funding

This work was supported by the National Institutes of Health (NIH R24 DK082841 to E.L.F. and J.H.); the American Diabetes Association (to E.L.F.); the Novo Nordisk Foundation (to E.L.F.); the Neuronetwork for Emerging Therapies at the University of Michigan, the A. Alfred Taubman Medical Research Institute (to S.A.E., S.E.E., and E.L.F.); the Applied Systems Biology Core of the George M. O’Brien Michigan Kidney Translational Core Center (NIH P30 DK081943); the NIDDK Diabetic Complications Consortium Pilot Grant (DiaComp, www.diacomp.org; DK076169; Sub-award #25034-75 to J.H.); the University of North Dakota (UND) Epigenetics Center of Biological Research Excellence Pilot Grant (CoBRE; NIH P20 GM104360 to J.H.); UND Post-Doc Pilot Grant (to K.G.); the multidisciplinary postdoctoral training program in basic diabetes research (T32DK101357 to S.E.E.); and the Nathan and Rose Milstein Research Fund (to S.A.E).

Author information

Affiliations

Authors

Contributions

K.G., S.A.E., S.E.E., E.L.F., and J.H. designed the studies. K.G., S.A.E., and C.P. conducted the studies. K.G., S.A.E., S.E.E., E.L.F., and J.H. analyzed the data. K.G., S.A.E., S.E.E., E.L.F., and J.H. wrote the manuscript. The authors read and approved the final manuscript.

Corresponding author

Correspondence to Junguk Hur.

Ethics declarations

Ethics approval and consent to participate

The use of human sural nerve samples was approved by the University of Michigan Institutional Review Board for Human Patient Research.

Consent for publication

Not applicable.

Competing interests

The authors have nothing to disclose.

Additional information

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

The original version of this article was revised: we were notified of a mistake in the way the author corrections have been implemented: a few random numbers had been incorrectly assigned to the second column in Figure 5A. These have now been removed.

Supplementary information

Additional file 1: Figure S1.

Sample distribution based on transcriptomic data. Figure S2. MA plot showing the top 50 DEGs in Group 1 versus Group 2. Figure S3. Gene expression patterns of the most significant 100 DEGs between Group 1 and Group 2. Figure S4. Functional enrichment analysis of DEGs by Reactome and DO. Figure S5. Functional enrichment analysis of DMGs by Reactome and DO. Figure S6. Functional enrichment analysis of overlapping DEGs and DMGs in the same direction.

Additional file 2: Table S1.

Patient characteristics. Table S2. RNA-seq & RRBS read quality and mapping summary. Table S3. Differentially expressed genes and their annotation. Table S4. KEGG enrichment analysis of DEGs. Table S5. DEG GO Clustering. Table S6. Reactome enrichment analysis of DEGs. Table S7. Disease ontology enrichment analysis of DEGs. Table S8. Differentially methylated CpG sites and their annotation. Table S9. KEGG enrichment analysis of DMGs. Table S10. DMG GO Clustering. Table S11. Reactome enrichment analysis of DMGs. Table S12. Disease ontology enrichment analysis of DMGs. Table S13. DEG & DMG overlap genes information.

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated in a credit line to the data.

Reprints and Permissions

About this article

Verify currency and authenticity via CrossMark

Cite this article

Guo, K., Eid, S.A., Elzinga, S.E. et al. Genome-wide profiling of DNA methylation and gene expression identifies candidate genes for human diabetic neuropathy. Clin Epigenet 12, 123 (2020). https://doi.org/10.1186/s13148-020-00913-6

Download citation

Keywords

  • Type 2 diabetes
  • Diabetic peripheral neuropathy
  • Transcriptomics
  • Epigenetics
  • DNA methylation
  • Gene expression
  • HbA1c
  • Human