Genome-wide DNA methylation profiling of human diabetic peripheral neuropathy in subjects with type 2 diabetes mellitus

ABSTRACT

DNA methylation is an epigenetic mechanism important for the regulation of gene expression, which plays a vital role in the interaction between genetic and environmental factors. Aberrant epigenetic changes are implicated in the pathogenesis of diabetes and diabetic complications, but the role of DNA methylation in diabetic peripheral neuropathy (DPN) is not well understood. Therefore, our aim in this study was to explore the role of DNA methylation in the progression of DPN in type 2 diabetes. We compared genome-wide DNA methylation profiles of human sural nerve biopsies from subjects with stable or improving nerve fibre counts to biopsies from subjects with progressive loss of nerve fibres. Nerve fibre counts were determined by comparing myelinated nerve fibre densities between an initial and repeat biopsy separated by 52 weeks. Subjects with significant nerve regeneration (regenerators) and subjects with significant nerve degeneration (degenerators) represent the two extreme DPN phenotypes. Using reduced representation bisulfite sequencing, we identified 3,460 differentially methylated CpG dinucleotides between the two groups. The genes associated with differentially methylated CpGs were highly enriched in biological processes that have previously been implicated in DPN such as nervous system development, neuron development, and axon guidance, as well as glycerophospholipid metabolism and mitogen-activated protein kinase (MAPK) signalling. These findings are the first to provide a comprehensive analysis of DNA methylation profiling in human sural nerves of subjects with DPN and suggest that epigenetic regulation has an important role in the progression of this prevalent diabetic complication.

KEYWORDS:

• Type 2 diabetes
• diabetic peripheral neuropathy
• DNA methylation
• reduced representation bisulfite sequencing
• RRBS

Introduction

In 2017, it was estimated that 451 million people worldwide suffer from diabetes, including 30 million Americans [1,2], with type 2 diabetes (T2D) accounting for approximately 90% to 95% of the diagnosed cases. T2D is primarily characterized by hyperglycaemia, hyperlipidaemia, and insulin resistance [3]. Both subjects with type 1 (T1D) and T2D experience diabetes-related microvascular complications which have a negative impact on the quality of life and affect multiple tissues including the nerve (diabetic neuropathy), kidney (diabetic nephropathy) and eye (diabetic retinopathy) [4]. Diabetic neuropathy is the most common microvascular complication of both T1D and T2D and can present in multiple ways, including small-fibre predominant neuropathy or autonomic neuropathy. However, diabetic peripheral neuropathy (DPN), a distal symmetric small and large fibre sensory-predominant neuropathy, is the most common type of diabetic neuropathy affecting the peripheral nervous system [5,6]. Affecting an estimated 50% of the patients with diabetes, DPN is a chronic, symmetrical, progressive disorder of the peripheral nervous system with early symptoms of pain, allodynia, and paraesthesias [7]. DPN initially affects the longest axons in the extremities and progresses in a distal to proximal, ‘stocking-glove’ pattern. The principal pathological abnormalities in DPN are axonal degeneration and segmental demyelination [8]. The sural nerve is a sensory nerve in the calf and quantitation of myelinated fibre density (MFD) of the sural nerve is considered one of the most reliable criteria for neuropathy diagnosis [9].

The pathogenesis of T2D and DPN is largely driven by lifestyle factors such as limited physical activity and poor diet [10]. However, there is a clear polygenic risk to developing T2D and DPN [11,12], and growing evidence suggests that in response to behavioural and environmental stressors, epigenetic mechanisms also contribute significantly to diabetes and its complications [13–15]. Indeed, epigenetic modifications are hypothesized as a way to bridge the gap between lifestyle factors and genetics [16]. In contrast to genetic alterations, epigenetic mechanisms regulate gene expression without changes in the underlying DNA sequence [17]. Aberrant epigenetic modifications, including DNA methylation, post-translational modifications of histones, chromatin remodelling, and non-coding RNAs [18], result in the dysregulation of up to 10% of the genes in humans [18,19].

DNA methylation refers to the covalent attachment of a methyl group to a cytosine adjacent to a guanine (CpG) by DNA methyltransferases [20]. While most commonly occurring in the promoter region of a gene, DNA methylation can also occur in various other genetic locations, including exons, introns, or even non-transcribed intergenic regions [21]. DNA methylation may have different influences on gene expression dependent upon where methylation takes place in the gene sequence [22,23]. DNA methylation is involved in the regulation of gene expression, maintenance of DNA structure, and control of transposable elements [24]; therefore, it plays an important role in many biological functions, including tissue differentiation and disease susceptibility [25,26]. Differential DNA methylation is also implicated in T2D development, insulin resistance, and obesity and may contribute to an intermediate stage of T2D pathogenesis [18,27,28]. Altered DNA methylation patterns are observed in multiple tissues of subjects with T2D compared to healthy controls [22,29–31]. Additionally, genome-wide and gene-specific studies support the idea that aberrant DNA methylation contributes to macro- and microvascular complications in diabetic subjects [32–34], and changes in DNA methylation are reported in other diabetic microvascular complications, such as in the kidneys of subjects with diabetic nephropathy [33,35] and in the blood of patients with diabetic retinopathy [34]. Furthermore, DNA methylation is associated with DPN in the transcriptional analysis of mouse models [3]. However, genome-wide DNA methylation patterns in subjects with DPN remain unexplored.

In this study, our aim was to explore the role of DNA methylation in the progression of DPN in T2D. Using reduced representation bisulfite sequencing (RRBS), we compared genome-wide DNA methylation in sural nerve biopsies obtained from subjects with T2D and either stable or improving DPN to subjects with progressive and worsening DPN. Our study provides the first comprehensive analysis of nerve-specific DNA methylation profiles and identifies a role for epigenetic regulation of previously implicated genes and pathways in nerve degeneration in T2D and DPN. Additionally, it identifies possible targets for the development of much-needed therapies for this common, morbid diabetic complication.

Results

Study population and genome-wide DNA methylation profiling

To identify pathways that are altered by epigenetic regulation in DPN, we first established genome-wide methylation patterns using RRBS in sural nerves samples from 12 DPN subjects from a larger double-blind placebo-controlled clinical trial where there was no effect of treatment [36,37]. Sural nerve samples were from subjects who had significant nerve degeneration (degenerators, n = 6) or regeneration (regenerators, n = 6) over a 52-week period as measured by myelinated fibre density (MFD). Subject characteristics and neuropathy phenotyping are presented in Tables 1 and 2. In total, an average of 88.2 million aligned reads were generated for each sample, with read lengths up to 48 bps, of which about 51.5 million could be uniquely mapped to the hg19 reference genome (unique mapping rate ranging from 56% to 60.9%). The overall percentage of methylated CpG sites ranged between 32.2% and 41.0% across the 12 samples (Supplemental Table S1), with no significant difference between degenerators and regenerators.

Table 1. Subject characteristics.

Group	Degenerator	Regenerator	P-value
Age (years)	53 ± 7	57 ± 8	0.385
BMI (kg/m^2)	28.9 ± 3.5	29.3 ± 4.2	0.851
Diabetes duration (years)	5.9 ± 5.8	12.1 ± 7.2	0.129
HbA1c (%)	7.6 ± 1.6	8.2 ± 1.7	0.561
Triglyceride (mmol/L)	11.0 ± 21.0	2.9 ± 1.2	0.387
Cholesterol (mmol/L)	7.0 ± 4.3	5.9 ± 0.8	0.559

BMI: body mass index; HbA1c: haemoglobin A1c. All subjects were diagnosed with T2D and were male. P-values were calculated by T-Test.

Table 2. Subject neuropathy phenotyping.

Group	Degenerator	Regenerator	P-value
MFD – baseline (fibres/mm^2)	2,954.6 ± 1,936.5	1,738.9 ± 950.7	0.198
MFD – 52 weeks (fibres/mm^2)	1,347.9 ± 1,052.6	2,830.4 ± 1,377.7	0.063
MFD change (fibres/mm^2)	−1,606.8 ± 1,118.3	1,091.5 ± 469.7	0.000
MFD percent change (%)	−53.4 ± 13.0	65.6 ± 18.1	0.000
O’Brien neuropathy score – baseline	3048.4 ± 1385.3	2521.4 ± 852.1	0.446
O’Brien neuropathy score – 52 weeks	3086.0 ± 1245.0	2436.8 ± 578.5	0.274

MFD: myelinated fibre density. O’Brien neuropathy score: a composite mean sum score of nerve conduction attributes, where the score is inversely correlated with the severity of neuropathy [80]. P-values were calculated by T-Test.

Sample information, alignment, and final cytosine methylation summary

Sample read coverage distribution per base as determined by methylKit [38] analysis of sequence alignment/map (SAM) data confirmed that the results had little to no PCR duplication bias, as no secondary peak was observed near the right side of the histograms (Supplemental Figure S1). A total of 1,527,600 CpG sites (excluding CpG sites in chromosomes X and Y) were identified across all 12 samples. Methylated CpG sites were classified by methylation levels into five sets: <5%, <10%, 10–90%, >90%, and >95% CpGs (Supplemental Table S2). More than half of the methylated CpG sites in each sample showed less than 10% methylation, regardless of the subject groups (P = 0.9). However, there was a significant difference in CpG sites with >95% methylation between the two groups (P = 0.0085), where the degenerators had greater overall hyper-methylation of CpG sites compared to regenerators (Figure 1). The DNA methylation clustering patterns identified by Ward’s hierarchical clustering methods illustrate that regenerator samples, and to a lesser extent degenerator samples, were tightly clustered, suggesting a distinct methylation footprint in regenerators vs. degenerators (Supplemental Figure S2).

Figure 1. Distribution of percentage methylation within groups. This histogram illustrates the overall distribution of the percent methylation of the measured 1,527,600 CpG sites. The boxes within each panel represent zoomed-in images covering the percent methylation between 90% and 100%.

Differential CpG methylation in degenerator versus regenerator sural nerve tissue

Using methylKit analysis of SAM data, we next determined the locations of differential DNA methylation in sural nerves of subjects with DPN. We identified a total of 3,460 differentially methylated CpGs (DMCpGs; 1,913 hypo-methylated and 1,547 hyper-methylated) using the regenerator group as reference (Supplemental Table S3). The numbers of DMCpGs per chromosome ranged between 67 and 293; however, there was no difference in DMCpG distribution (per one million base pairs) between chromosomes (Supplemental Figure S3 and supplemental Table S4).

Effects on gene expression are dependent on the location of DNA methylation within the gene sequence [22]. We therefore next determined which parts of the genome were affected by altered methylation using genomation [39] tool. Approximately 21.2% of DMCpGs were located in promoter regions, 9.9% in exons, 36.1% in introns, and 32.8% in intergenic regions (Figure 2(a)). The distribution of the DMCpG sites in relation to CpG islands (regions along the genome with a high number of CpG sites; GC fractions >0.5, lengths >200 bp, and a ratio of observed-to-expected CpG > 0.6) is shown in Figure 2(b). We found 130 CpG sites (3.8%) in CpG islands, and 331 CpG (9.6%) sites in CpG shores (regions within 2 kb of a CpG island). Moreover, approximately 86.7% of DMCpGs were found in non-CpG-rich regions, which were >2 kb from a CpG island (denoted as ‘other’ in the figure).

Figure 2. Distribution of DMCpGs (degenerators compared to regenerators) according to the gene region. (a) Pie chart showing the percentage of DMCpGs located in promoters, exons, introns, and intergenic regions. (b) Pie chart showing the percentage of DMCpGs located in CpG islands (CpGi), CpG shores (regions within 2 kb of a CpG island), and others.

Annotating differential methylation events

To better understand how DMCpGs may result in functional changes in the genome, using genomation analysis DMCpGs were mapped to 2,835 unique genes on the basis of their genomic location to identify hypo- or hyper-methylated genes (Supplemental Table S3). Among the observed DMCpGs, 1,913 corresponded to 1,607 unique genes exhibiting hypo-methylation, while 1,547 DMCpGs corresponded to 1,352 unique genes exhibiting hyper-methylation in degenerator samples compared to regenerator samples. The top 15 most differentially methylated CpGs and their annotated genes (DMGs) are presented in Table 3 and include multiple non-coding and microRNAs (miRNAs) as well as coding genes. Among all differentially methylated genes (DMGs), 102 genes (represented by 102 Ensemble transcript IDs) did not have valid official gene symbols and were therefore removed from further analyses. Of the remaining genes, 516 included two or more DMCpGs. Moreover, 1,295 out of 2,835 genes had DMCpGs located within a 5 kb distance from their transcript start sites (TSSs) (DMCpGs located primarily within promotor and enhancer regions).

Table 3. Top 15 differentially methylated CpGs and their annotated genes.

Gene	CpG location	Q-Value	Methylation difference	Type	Description
RNU7-88P	Chr13:59439918	1.60E-116	−68.47	snRNA	RNA, U7 small nuclear 88 pseudogene
CFAP46	Chr10:132914032	8.23E-109	62.75	Protein Coding	cilia and flagella associated protein 46
AC114814.4	Chr2:235044515	1.87E-101	62.40	LincRNA
MIR3138	Chr4:10079577	6.57E-97	60.31	microRNA	miR3138
CASZ1	Chr1:10768122	1.25E-89	54.57	Protein Coding	castor zinc finger 1
SLC38A10	Chr17:81287547	3.64E-88	−47.39	Protein Coding	solute carrier family 38 member 10
ACPP	Chr3:132302985	3.07E-74	62.19	Protein Coding	acid phosphatase, prostate
PTGIS	Chr20:49560881	9.13E-74	−52.41	Protein Coding	prostaglandin I2 (prostacyclin) synthase
DHX37	Chr12:124970826	2.46E-70	48.09	Protein Coding	DEAH-box helicase 37
SHANK2	Chr11:70460640	2.37E-68	53.75	Protein Coding	SH3 & mult. ankyrin repeat domains 2
MIR885	Chr3:10379191	5.57E-68	−54.63	microRNA	miR885
RP11-302L19.1	Chr6:170134734	6.38E-68	−61.26	ncRNA
AC009531.2	Chr7:35363380	1.03E-67	−59.35	LincRNA
KDM4D	Chr11:94998595	1.49E-67	−53.79	Protein Coding	lysine demethylase 4D
ASPG	Chr14:104101239	1.00E-66	35.84	Protein Coding	asparaginase

CpG location corresponds to the coordinate of human genome hg19 assembly. Methylation difference: the difference of methylation (%) between the two groups.

Enriched biological functions among DMCpGs

To determine potential biological functions of the identified DMCpGs and their corresponding DMGs, enrichment analysis was performed on DMGs using Gene Ontology (GO) terms and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways. Using GO analysis, DMGs in degenerators compared to regenerators were primarily involved in biological processes related to the nervous system, such as ‘nervous system development’, ‘neuron development’, and ‘neuron differentiation’ (Figure 3; Supplemental Table S5). Within these pathways, top DMGs included genes highly involved in nerve function and repair, such as netrin-4 and dihydropyrimidinase like 2.

Figure 3. GO enrichment analysis of sural nerve DMGs (degenerators compared to regenerators). The top 20 most significantly enriched GO terms among DMGs are shown. Rich factor is the proportion of DMGs to all genes belonging in each GO term. The values on each bar correspond to the numbers of DMGs in each GO term, and colour gradient (see colour scale right of the figure) indicates the level of significance.

KEGG analysis showed a significant enrichment of pathways know to be involved in diabetes and DPN, such as ‘MAPK signalling’ and ‘glycerophospholipid metabolism’, as well as nervous system pathways including ‘axon guidance’ (Figure 4; Supplemental Table S6). Significant DMGs that were highly different between degenerators and regenerators within these pathways included genes related to lipid metabolism and diabetes, such as: phospholipase A2, phosphatidylserine decarboxylase, and vascular endothelial growth factor.

Figure 4. KEGG pathway enrichment analysis of sural nerve DMGs (degenerators compared to regenerators). All significantly enriched KEGG terms among DMGs are shown. Rich factor is the proportion of DMGs to all genes belonging in each KEGG pathway. The size of dot corresponds to the numbers of DMGs in each KEGG pathway, and colour gradient (see colour scale right of the figure) indicates the level of significance.

Validation of DMCpGs and DMGs

To perform a targeted validation of the results obtained by the genome-wide DNA methylation analysis, we quantified DNA methylation by pyrosequencing. Three CpG sites (Supplemental Table S7) located in promoter regions (< 1kb from TSSs) were selected based on functional enrichment analysis (Supplemental Tables S5 and S6) or degree of differential methylation from a list of the top 30 DMGs (Supplemental Table 3). Of the three selected, dihydropyrimidinase like 2 (DPYSL2; a gene implicated in nervous system development) showed a significantly higher degree of methylation in the degenerator group compared to the regenerator group by pyrosequencing (p = 0.0005), and similar results were observed for miR3138, the most differentially methylated miRNA in our analysis (p = 0.0130) (Figure 5(a)). For DPYSL2 and miR3138, the fold difference between the regenerator mean percent methylation and the degenerator samples presented a similar pattern between genome-wide RRBS and loci-specific pyrosequencing methylation analyses. In addition, their gene expression as measured by RT-PCR followed a decreasing trend compared to their increased methylation (Figure 5(b)). In contrast, significant changes in methylation between the two groups for the nuclear factor of activated T cells 1 (NFATC1) were observed only with RRBS but not with pyrosequencing (Figure 5(a)).

Figure 5. Cross-validation of methylation profiling and sequence-specific methylation of DYSPL2 and miR3138. (a) Percent CpG methylation for DPYSL2, miR3138, and NFATC1 measured by conventional pyrosequencing. The data are presented as mean ± SEM, and an unpaired t-test was used. (b) Methylation and gene expression fold change in the degenerators relative to the regenerators within analysis platform (RRBS, pyrosequencing, and qPCR). Dashed line represents the average gene expression by qPCR analysis of the regenerator group. *p < 0.05, **p < 0.01, ***p < 0.001. For RRBS: DEGEN n = 6; REGEN n = 4. For pyrosequencing: DEGEN n = 5; REGEN n = 4. For qPCR: DEGEN n = 4, REGEN n = 3.

Discussion

We have previously reported changes in gene expression in nerves from mouse models [3,40,41] as well as from human subjects with DPN [42]. In the current study, we examined the epigenetic modifications of human sural nerves from subjects with varying degrees of DPN. RRBS was used to determine differential methylation in DPN subjects who had significant sural nerve degeneration over a one-year period compared to subjects who had stable disease or significant regeneration [4]. Using this technique, we identified 3,460 DMCpGs mapped to 2,835 unique genes. Among the DMCpGs, 1,913 (with 1,607 unique genes) were hypo-methylated and 1,547 (with 1,352 unique genes) were hyper-methylated. Functional analysis indicated that epigenes were enriched in biological pathways related to nervous system development, axon guidance, glycerophospholipid metabolism, and MAPK signalling; several of which have previously been shown to be involved in DPN [42–44]. Further, we validated our results using pyrosequencing and RT-PCR, confirming changes in methylation as well as a trend for a change in gene expression of DPYSL2 and miR3138. Our results suggest that gene expression and associated pathophysiological pathways implicated in DPN progression may be under epigenetic control.

DNA methylation occurs in various regions along the gene [45]. The sites of DNA methylation can effect gene transcription, for example, gene transcription can be halted by methylation within promoter regions and activity can be increased when methylation occurs within intragenic regions [23]. However, while aberrant promoter methylation is a hallmark of diseases such as cancer and is associated with gene silencing [46], changes in DNA methylation commonly occur in regulatory regions rather than in promoter regions in diabetes and diabetic kidney disease [35,47]. These findings are consistent with the high percentage of DMCpGs in intergenic and intronic regions observed in our study. Changes in DNA methylation within regulatory regions may also affect the expression of long non-coding RNAs or miRNAs. MiRNAs themselves play an important role in health and disease [48]. Moreover, epigenetic regulation of miRNAs, including regulation via methylation, is involved in neurodegenerative disorders [49]. In the current study, we observed multiple DMCpGs mapped to annotated miRNAs and long non-coding RNAs.

To better understand the potential biological implications of changes in DNA methylation with progressive DPN, we performed pathway enrichment of significant DMCpGs using GO and KEGG analysis. Given that subject groups were stratified based on sural nerve degeneration vs. regeneration, it is not surprising that pathway analysis revealed multiple significantly enriched pathways related to nervous system development and axon guidance. Within these enriched pathways, a significant DMG in the current study was netrin-4 (NTN4), a gene which promotes neurite elongation and outgrowth [50,51]. In a report comparing methylation profiles in subjects with diabetic nephropathy, differential methylation of saliva netrin-4 was associated with end-stage renal disease [33]. In the nervous system, netrin-4 plays a key role in thermal and mechanical sensitivity during neuropathic pain [52], and we have previously reported an association between peripheral nerve Ntn4 gene expression and DPN phenotype in the db/db mouse model of T2D and DPN [40,41,44]. Taken together, our data suggest that differential DNA methylation of NTN4 regulates its gene expression, which in turn plays an important role in nerve regeneration and repair in DPN.

Glycerophospholipid metabolism was another significantly enriched pathway identified between degenerators and regenerators. Glycerophospholipids are major components of neural membranes and changes in their levels are linked to neurological disorders [53]. Our findings suggest that differential DNA methylation of glycerophospholipid metabolism pathway genes contributes to DPN pathogenesis and disease progression. Importantly, particular DMGs of interest within the glycerophospholipid metabolism pathway were those previously associated with nerve injury, including phospholipase A2 and phosphatidylserine decarboxylase. Phospholipase A2 is involved in inflammation-mediated metabolic disorders, including diabetes and obesity [54]. Phospholipase A2 is also required for neurite growth in the peripheral nervous system [55] and is involved in myelin breakdown during Wallerian degeneration [56]. Phosphatidylserine decarboxylase, an enzyme that catalyses the reaction from phosphatidylserine to phosphatidylethanolamine, is associated with neurodegeneration [57,58]. Further, phosphatidylethanolamine supplementation prevents myelination deficits and restores nerve function in a Charcot-Marie-Tooth disease 1A animal model [59]. Collectively, our results coupled with the established roles of glycerophospholipids in nerve function highlight important roles for phospholipase A2 and phosphatidylserine decarboxylase that require further investigation [57,60].

Another significantly enriched pathway in the degenerator group compared to regenerators was ‘MAPK signalling’. The MAPK signalling pathway responds to a variety of signals and has roles in regulating cell death, differentiation, growth, and inflammation [61]. In addition, it modulates central and peripheral nociceptive information in animal models of neuropathic pain [62], and is well recognized in many neurodegenerative disorders [63], including DPN [64]. Furthermore, the three arms of the MAPK signalling family, ERK (extracellular signal-regulated kinases), JNK (c-Jun N-terminal kinases), and p38 are altered in multiple DPN models as well as in humans with DPN [65,66]. Additionally, ERK, JNK, and p38 are critical for myelination of the axons and are involved in nervous system degeneration and regeneration [64]. In the current study, the enrichment of MAPK DMGs in the regenerator group supports a role for MAPK signalling in DPN progression and suggests that MAPK enzymes may be potential disease-modifying therapeutic targets.

We performed pyrosequencing and RT-PCR validation of DPYSL2 (Dihydropyrimidinase Like 2; also known as CRMP2), NFATC1 (Nuclear Factor Of Activated T Cells 1), and miR3138. DPYSL2 was one of the genes chosen for pyrosequencing as it was highly significantly different between groups within both GO and KEGG pathway analyses. Studies have implicated DPYSL2 in the regulation of neuron guidance, growth, polarity [67], and synaptic signalling [68,69]. Additionally, altered DPYSL2 levels are associated with neurodegenerative diseases, as well as with neuropathic pain [69–71]. We previously reported a significant decrease in the gene expression of Dpysl2 in the sciatic nerves of db/db neuropathic mice compared to normal controls [43]. As expected, differential methylation of DPYSL2 was confirmed with pyrosequencing and was associated with a change in gene expression as measured by RT-PCR. Our results indicate that differential DNA methylation is a potential regulator of DPYSL2, a gene with a known role in nervous system function and neuropathy.

Also chosen for validation was the most highly differentially methylated miRNA between the degenerators and regenerators, miR3138, a miRNA involved in tumour suppression [72]. While little is known about its role in diabetes and diabetic complications, its predicted targets (according to the target prediction database, miRDB [73]) include DMGs or gene families observed in our results, including the immune system related nuclear factor of activated T-cells and neurexin 2, a gene encoding for cell adhesion molecules. Additionally, we cross-referenced predicted and known targets of this miR3138 from miRDB and DIANA-TarBase (database of experimentally validated microRNA target) [74] against the differentially expressed genes between progressive (similar to degenerator) and non-progressive (similar to regenerator) DPN from our previous microarray-based gene expression study in human sural nerves [42]. Using this approach, we identified erb-b2 receptor tyrosine kinase 4 (ERBB4) as a target of miR3138 and a significant gene involved in the progression of DPN [42]. ERBB4 has also previously been associated with sensory nerve loss and nerve regeneration [75,76], providing evidence for a potential regulatory role of miR3138 in nerve degeneration and regeneration in DPN. We confirmed that this differential methylation seen by RRBS and pyrosequencing was associated with a trending decrease in miR3138 gene expression in degenerator sural nerves. These data suggest that epigenetic downregulation of miR3138 expression may play a role in DPN progression. Further studies are needed to fully understand its role in DPN.

Our results suggest that DNA methylation is altered in sural nerves from subjects with DPN, and that the state of methylation may contribute to the regenerative ability of the nerve. Limitations of this study include sample size, and the limited scope of RRBS (as it interrogates only about 10% CpGs of the human genome [77]). Additionally, while two of the three chosen sites were confirmed using pyrosequencing, one (NFATC1) failed. This may be due to the number of CpG sites, as more heavily populated CpG regions did not pass pyrosequencing quality control or validation. In particular, the amplicon analysed for NFATC1 contained up to six potential CpG sites, and homopolymeric sequences limit accurate sequence resolution with genome-wide methylation profiling [78], and therefore may have introduced ambiguity in the RRBS measurements. Lastly, the current analysis did not account for cell-type heterogeneity of sural nerves that may be altered between the two study groups. Adjusting for changes in cell-type composition would provide new insights into the cell-specific processes that drive neuropathy pathogenesis. Future experiments are warranted to investigate functional and cell-specific changes due to differences in methylation in highly methylated regions and connection to gene expression at corresponding loci in a larger cohort of patients.

In summary, our study is the first to examine changes in DNA methylation in human sural nerves of subjects with DPN, and links specific DNA methylation profiles with DPN pathogenesis, specifically disease progression. In particular, DNA methylation contributes to the regulation of functional pathways, such as nervous system development, axon guidance, glycerophospholipid metabolism, and MAPK signalling. As such, it provides novel and important insights into epigenetic regulation of critical disease-modifying pathways that can be targeted to understand DPN pathogenesis and develop mechanism-based therapies.

Materials and methods

Subject selection and characteristics

Biopsies from human sural nerves were obtained during a previous 52-week double-blind placebo-controlled clinical trial of acetyl-L-carnitine treatment for DPN [36,37]. Criteria for enrolment in the larger cohort from which the samples for this manuscript were selected is detailed in a previously published study [36]. Briefly, adult male and female subjects diagnosed with T1D or T2D for a minimum of 1 yr and a HbA1c > 5.9% were included. Additionally, all subjects were phenotyped for neuropathy by a minimum of two neurological findings, including clinical symptoms or abnormal electrophysiological tests (nerve conduction velocity or vibration perception threshold) [79,80]. Subjects with complicating diseases (such as HIV or significant cardiac or hepatic disorders), alcohol or drug abuse, or other causes of neuropathy were excluded. MFD of the sural nerve is considered the most reliable morphological criteria for DPN diagnosis [9]. Within this cohort, acetyl-L-carnitine treatment did not affect MFD, and we previously reported that there were separate groups of subjects with DPN who demonstrated either significant peripheral nerve degeneration or significant regeneration (denoted as degenerators and regenerators, respectively), as measured by MFD percent changes over 52 weeks [4]. Degenerators (n = 6) and regenerators (n = 6) with the largest changes in MFD were selected for the current methylation profiling study. All selected subjects were diagnosed with T2D and were male; subject characteristics are summarized in Table 1, and neuropathy phenotyping is presented in Table 2 (See Supplemental Tables S8 and S9 for the individual subject-level data). The use of human sural nerve samples was approved by the University of Michigan Institutional Review Board for Human Subject Research.

DNA extraction and RRBS

DNA was extracted from sural nerve tissue and prepared for RRBS by the University of Michigan Epigenomics Core. Briefly, DNA extraction followed the protocol published by Garrett-Bakelman et al. [81]. Sample quantity was measured using a Qubit fluorometer, and the quality of the genomic DNA (gDNA) was evaluated using the gDNA assay on Agilent’s TapeStation. High-quality gDNA was digested with the MspI restriction enzyme and purified using phenol:chloroform extraction and ethanol precipitation, before blunt-end digestion, phosphorylation, and ligation of an adapter duplex with methylated cytosines. The ligated fragments were cleaned and processed for size selection on an agarose gel. Selected fragments were bisulfite converted, followed by PCR amplification and cleanup with AMPure XP beads. The libraries were quantified using the Qubit assay and Agilent’s High Sensitivity D1000 assay, and then sequenced on the Illumina Hi-Seq 4000 platform to generate raw sequence RRBS reads.

DNA sequence processing and alignment

The raw sequence reads were evaluated and processed using software and methods as outlined previously [82]. Briefly, the reads were first assessed for quality using FastQC (version 0.11.5, Babraham Bioinformatics, UK), and then low-quality reads were trimmed off using Trim Galore (version 0.4.2, Babraham Bioinformatics, UK). Only the bases with a Phred > 30 were kept for downstream analysis and Illumina adaptors were also removed. A special setting designed for trimming RRBS data (–rrbs) in Trim Galore was also used. Specifically, Trim Galore trimmed the first two bases from the 3′ end of the reads so the cytosine (C) base closest to the second enzyme-cut site was not included in methylation calling, since the RRBS approach introduces artificial CpG sites at the end. Trimmed reads were then aligned and mapped with Bismark version 0.16.3 [83]/Bowtie2 version 2.2.9 [84] to an in silico bisulfite converted human reference genome (GRCh38). The unique alignment was then identified and used to make a methylation call after discarding reads with multiple mapping.

Methylome characterization and differential methylation calculation

Generated sequence alignment/map (SAM) format files were sorted and then directly read in R using the methylKit package version 1.2.4 [38]. Read coverage per base and percent methylation per base as the basic information were calculated and displayed, and the methylation level at each cytosine was estimated by #C/(#C+#T), where #C is the number of methylated reads and #T is the number of unmethylated reads. CpG sites with minimum read coverage of >10 and minimum Phred quality scores >20 were selected for further analysis. Genome-wide methylation patterns were compared using Pearson’s correlation test. Clustering analysis was carried out using Ward’s minimum variance method to determine methylation levels and patterns [85].

DMCpGs were defined as a > 25% difference in cytosine methylation level between regenerators and degenerators, with a q-value from MethylKit, equivalent to false discovery rate (FDR) adjusted p-value, <0.01. Briefly, p-values of DMCpGs were calculated using logistic regression by comparing the number of methylated reads (Cs) and unmethylated reads (Ts) between degenerators and regenerators, and then corrected for multiple testing using FDR. For DMR analysis, a sliding windows approach with 1 kb width was used to summarize the methylation information on those tiles. Then, the same significance criteria were used to determine DMRs (>25% change with a q-value < 0.01).

Filtered DMCpGs and DMRs were annotated based on genes and CpG island features, including gene bodies and 2, 5, and 10 kb regions upstream from TSSs using genomation [39] R package. CpG islands were identified by three criteria; GC fractions > 0.5, lengths > 200 bp, and a ratio of observed-to-expected CpG > 0.6. CpG shores were identified by length (2 kb) and adjacent positions to CpG islands. Annotation of human CpG islands was obtained from the University of California, Santa Cruz (http://genome.ucsc.edu/; hg38), and annotation of the human genes was obtained from Ensembl (http://www.ensembl.org/; protein-coding genes).

Functional enrichment analysis

To identify common biological functions among the nearest genes of DMCpGs, deemed as DMGs, and their promotor regions (1 kb upstream of TSSs), gene set enrichment analysis was performed to identify significantly over-represented biological functions in terms of GO (http://www.geneontology.org) biological process terms and KEGG (http://www.genome.jp/kegg/) pathways. Hypergeometric test was used to examine the statistical significance of over-representation, with an adjusted p-value of < 0.05 as cutoff.

Pyrosequencing

Design of primers for DPYSL2, NFATC1, and miR3138 and preparation of the pyrosequencing amplicons were done at the University of Michigan Epigenomics Core using a published protocol [86]. Briefly, 20 ng of bisulfite-treated DNA was used for PCR amplification (45 cycles) using HotStarTaq Master Mix Kit (203443; Qiagen) and 200 pM of amplicon-specific primers (Supplemental Table S7). For DYSPL2, 5% DMSO was added to the PCR reaction. To verify the expected size of the amplicons, 5 μl of each PCR reaction was run on an ethidium bromide-agarose gel. The amplicons were submitted to the University of Michigan Sequencing Core for sample preparation and pyrosequencing on the PyroMark Q96 MD [87].

Gene expression

Sural nerve RNA was isolated using a RNeasy fibrous tissue mini kit (74704; Qiagen). cDNA template was generated using 200 ng of total RNA and iScript Supermix (1708840; BioRad) in a 40 μl reaction. The reaction was diluted 1:1 with dH₂O. qPCR reactions were performed in triplicate using TaqMan™ 2X gene expression Master Mix (4369016; ThermoFisher/Applied Biosystems), 2 μl of template, and sequence-specific TaqMan™ probes (ThermoFisher/Applied Biosystems; Supplemental Table S7) in an Applied Biosystems StepOne^TM RT-PCR system. The conditions used were: 40x cycles of 30 s at 95º C, 60 s at 55–60º C, and 30 s at 72º C. This was followed by a 5 min final phase at 72º C. C_T values were used to calculate ΔC_T and ΔΔC_T (with YHWAZ as a housekeeping control and regenerators as the relative control group; Supplemental Table S7). Data are expressed as the mean of the relative quantity of gene expression (2^−ΔΔCT) [88]. Statistically significant differences were calculated using a t-test with GraphPad Prism 7 software (GraphPad software Inc.).

Acknowledgments

The authors thank Dr. Benjamin Murdock and Dr. Stacey Sakowski-Jacoby at the University of Michigan for their expert editorial advice. Some of the data were presented as an abstract at the Keystone Bioenergetics and Metabolic Disease Symposia in Keystone, CO in January 2018.

Supplementary material

Supplemental data for this article can be accessed here.

Additional information

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 Program for Neurology Research and Discovery at University of Michigan, the A. Alfred Taubman Medical Research Institute, 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 [CoBRE; NIH P20 GM104360]; Pilot Grant (to J.H.), UND Post-Doc Pilot Grant (to K.G.), the multidisciplinary postdoctoral training program in basic diabetes research [T32DK101357; to S.E]; and the Nathan and Rose Milstein Research Fund.