Genome-Wide DNA Methylation and LncRNA-Associated DNA Methylation in Metformin-Treated and -Untreated Diabetes

Metformin, which is used as a first line treatment for type 2 diabetes mellitus (T2DM), has been shown to affect epigenetic patterns. In this study, we investigated the DNA methylation and potential lncRNA modifications in metformin-treated and newly diagnosed adults with T2DM. Genome-wide DNA methylation and lncRNA analysis were performed from the peripheral blood of 12 screen-detected and 12 metformin-treated T2DM individuals followed by gene ontology (GO) and KEGG pathway analysis. Differentially methylated regions (DMRs) observed showed 22 hypermethylated and 11 hypomethylated DMRs between individuals on metformin compared to screen-detected subjects. Amongst the hypomethylated DMR regions were the SLC gene family, specifically, SLC25A35 and SLC28A1. Fifty-seven lncRNA-associated DNA methylation regions included the mitochondrial ATP synthase-coupling factor 6 (ATP5J). Functional gene mapping and pathway analysis identified regions in the axon initial segment (AIS), node of Ranvier, cell periphery, cleavage furrow, cell surface furrow, and stress fiber. In conclusion, our study has identified a number of DMRs and lncRNA-associated DNA methylation regions in metformin-treated T2DM that are potential targets for therapeutic monitoring in patients with diabetes.


Introduction
DNA methylation, the most widely studied epigenetic mechanism, involves the covalent addition of a methyl group at the 5 position of the cytosine ring within the 5 -CpG-3 dinucleotides to create a 5-methylcytosine (5-mC). The target of DNA methylation, catalyzed by DNA methyltransferases (DNMTs) enzymes, are CpG nucleotides, which are usually unmethylated [1]. These CpG nucleotides occur at high-frequency in the promoter regions of genes and are frequently associated with hyper-or hypomethylation events [2]. Hypermethylation of promoter CpG islands can result in suppression of gene expression, whereas hypomethylation is associated with the transcriptional activation of affected genes [3]. Various studies suggest that these modifications may alter the transcriptional activity of genes and contribute to pathogenic conditions, such as the type 2 diabetes mellitus (T2DM) A total of 33 differentially methylated regions (DMRs) were observed between individuals on metformin treatment compared to screen-detected subjects. Of these, 22 were hypermethylated, whilst another 11 were hypomethylated in participants treated with metformin, and these are summarized in Table 2. Lnc-associated DNA methylation peaks in the promoter regions are summarized in Table 3 showing that 36 were hypermethylated, and 21 were hypomethylated in individuals on metformin. KEGG pathway analysis revealed no enriched pathways. Based on GO analyses, we retrieved the biological process, cellular process, and molecular function of the DMRs, and these are presented in Figures 1 and 2. The top enrichment scores for cellular processes of hypermethylated DMRs in subjects on metformin were associated with the axon initial segment, node of Ranvier, cell periphery, cleavage furrow, cell surface furrow, and stress fiber (Figure 1), whilst the hypomethylated biological processes were associated with photoreceptor outer segment ( Figure 2).   Gene name refers to the name of the DMR-associated gene. Genomic coordinates refers to the genomic locus of the DMR. DMR Length refers to the length of the DMR. log2FC refers to the fold change of normalized tag counts between two groups (log2 transformed). The p-value refers to the p-value of the DMR, the smaller, the more significant. The q-value refers to BH FDR corrected p-value.

Discussion
In this study, we measured DNA methylation in diabetic individuals on metformin treatment compared to newly diagnosed diabetes cases and found 33 differentially methylated regions (DMRs) of which 22 (67%) were hypermethylated in diabetes subjects on metformin therapy. Furthermore, 57 lncRNA-associated DNA Methylation regions (36 hypermethylated and 21 hypomethylated) were detected of which 63% were hypermethylated in metformin-treated subjects. Functional pathway analysis of these DMRs revealed that they affect gene expression in the axon initial segment (AIS), node of Ranvier, cell periphery, cleavage furrow, cell surface furrow, and stress fiber.
Amongst the hypomethylated DMRs found in this study were genes in the SLC family, specifically SLC25A35 and SLC28A1. The SLC family is known for its importance in drug development, and their proteins include passive transporters, symporters, and antiporters and are located in cellular and organelle membranes [23]. Transporters facilitate the movement of a specific substrate across the membrane with or against its concentration gradient and sequence analysis of SLC25A35 indicates that it likely functions as an oxaloacetate carrier, implying mitochondrial association [24]. On the other hand, SLC28A1, a high-affinity pyrimidine nucleoside transporter, plays a role in renal reabsorption and has been observed to be impaired during diabetes [25]. Metformin treatment has been associated with lower methylation levels in SLC transporter genes, as was shown in a study conducted on metformin transporter genes in liver tissue [18]. Mitochondrial dysfunction due to diabetes affects oxidative phosphorylation and decreases ATP production. As SLC proteins transport various solutes across the mitochondrial membrane in order to partake in a number of metabolic pathways [26,27], the decrease in methylation and subsequent increase in gene expression of SLC transporters could be indicative of the antidiabetic effect of metformin treatment.

Discussion
In this study, we measured DNA methylation in diabetic individuals on metformin treatment compared to newly diagnosed diabetes cases and found 33 differentially methylated regions (DMRs) of which 22 (67%) were hypermethylated in diabetes subjects on metformin therapy. Furthermore, 57 lncRNA-associated DNA Methylation regions (36 hypermethylated and 21 hypomethylated) were detected of which 63% were hypermethylated in metformin-treated subjects. Functional pathway analysis of these DMRs revealed that they affect gene expression in the axon initial segment (AIS), node of Ranvier, cell periphery, cleavage furrow, cell surface furrow, and stress fiber.
Amongst the hypomethylated DMRs found in this study were genes in the SLC family, specifically SLC25A35 and SLC28A1. The SLC family is known for its importance in drug development, and their proteins include passive transporters, symporters, and antiporters and are located in cellular and organelle membranes [23]. Transporters facilitate the movement of a specific substrate across the membrane with or against its concentration gradient and sequence analysis of SLC25A35 indicates that it likely functions as an oxaloacetate carrier, implying mitochondrial association [24]. On the other hand, SLC28A1, a high-affinity pyrimidine nucleoside transporter, plays a role in renal reabsorption and has been observed to be impaired during diabetes [25]. Metformin treatment has been associated with lower methylation levels in SLC transporter genes, as was shown in a study conducted on metformin transporter genes in liver tissue [18]. Mitochondrial dysfunction due to diabetes affects oxidative phosphorylation and decreases ATP production. As SLC proteins transport various solutes across the mitochondrial membrane in order to partake in a number of metabolic pathways [26,27], the decrease in methylation and subsequent increase in gene expression of SLC transporters could be indicative of the antidiabetic effect of metformin treatment. It is, therefore, likely that metformin in its demethylation action of SLC mitochondrial carriers could possibly aid cell repair in these patients, however, this requires further investigation.
Functional pathway analysis observed in this study is consistent with the basic pathological abnormalities in Diabetic Peripheral Neuropathy (DPN), such as axonal degeneration and demyelination, lack of sensation, numbness, paresthesia, and allodynia experienced by diabetic individuals [28]. Cell death of nerves in DPN results from multifactorial metabolic imbalances associated with diabetes. The resulting mitochondrial dysfunction through a series of cascade effects involving AMP-activated protein kinase (AMPK), sirtuin (SIRT), and peroxisome proliferator-activated receptor-γ coactivator α (PGCα) suppresses mitochondrial oxidative phosphorylation, resulting in neuronal and axonal degeneration through increased oxidative injury [29,30].
Treatment with metformin was shown to decrease the incidence of DPN as was observed by the Bypass Angioplasty Revascularization Investigation 2 Diabetes trial [31]. Although metformin cannot reverse the nerve damage caused by diabetes, it could assist in managing blood glucose levels and improving the symptoms for patients.
In addition to DMRs, 57 lncRNA-associated DNA Methylation Peaks were detected when comparing known diabetic individuals to screen detected patients. Most recently the NONCODE database has updated the numbers of human lncRNAs to 167,150 with numbers still increasing [32]. Recent genome-wide association studies (GWAS) have shown positive correlation of some lncRNAs and diabetes [33]. In a related study, Sathishkumar et al. (2018) found increased levels of lncRNAs in T2DM patients, including HOTAIR, MEG3, LET, MALAT1, MIAT, CDKN2BAS1/ANRIL, XIST, PANDA, GAS5, Linc-p21, ENST00000550337.1, PLUTO, NBR2THRIL, and SALRNA1. The majority of these lncRNAs were involved in cell cycle regulation and senescence with their expression levels correlating to poor glycemic control, insulin resistance, and inflammation [11]. Similarly, HECTD4 and MBTPS1 were identified as the target genes for lncRNAs ENST00000364558 and ENST00000565382, respectively, with involvement in the development of T2DM by means of the lysosome and phagocytic signaling pathways [34]. Our findings indicate several novel lncRNA, including a lncRNA associated with the mitochondrial ATP synthase-coupling factor 6 (ATP5J) enzyme thought to be involved in the oxidative phosphorylation pathway [35]. Our data suggest higher methylation levels of this lncRNA in metformin-treated subjects, possibly pointing to suppression of this lncRNA allowing for ATP5J expression. Although little association was found between metformin and lncRNAs in our study, the significant novel lncRNA identified warrants further investigation to explore possible roles in type 2 diabetes.
The limitations of this study include the small sample size and the inclusion of women only; however, this allowed comparison and limited error that may result in statistical manipulation of a small sample size by sex. Furthermore, we used peripheral blood DNA to perform the genome-wide DNA methylation analysis. Epigenetic changes are believed to be organ specific; however, investigations on peripheral blood DNA have shown consistent methylation patterns with other organs [36,37]. Although the average (5.2 years) duration of disease in metformin-treated subjects was within the four to six years in which a person may have had the condition before clinical diagnosis [38], these findings should be interpreted with caution. In conclusion, our study has identified a number of DMRs and lncRNA-associated DNA methylation regions in metformin-treated T2DM that are potential targets for therapeutic monitoring in diabetes patients. However, these findings require further longitudinal study investigations that can clearly ascertain that these observations are not confounded by the duration and severity of diabetes.

Ethical Approval of the Study
This investigation used data from the Cape Town Vascular and Metabolic Health (VMH) study), which were approved by the Research Ethics Committees of the Cape Peninsula University of Technology and Stellenbosch University (resp., NHREC: REC-230 408-014 and N14/01/003; approved date: 21 May 2018). The Code of Ethics of the World Medical Association (Declaration of Helsinki) was also applied to the study. Signed written consent was obtained from all participants after all procedures were explained in the language of their choice.

Study Procedures
In this case-control study, the participants were females matched for both age and body mass index. All study participants underwent a standardized interview, blood pressure, and anthropometric measurements. A 75 g oral glucose tolerance test (OGTT) was performed on participants with no previous diagnosis of diabetes mellitus. Participants who met the World Health Organisation (WHO) criteria for diabetes were termed as screen-detected or newly diagnosed diabetes. Biochemical parameters analyzed at an ISO 15189 accredited pathology practice (PathCare, Reference Laboratory, Cape Town, South Africa) included the following: plasma glucose, serum insulin, serum creatinine, total cholesterol (TC), high-density lipoprotein cholesterol (HDL-c), triglycerides (TG), low-density lipoprotein cholesterol (LDL), C-reactive protein (CRP), γ-glutamyl transferase (GGT), AST, ALT, and glycated hemoglobin (HbA1c), certified by the National Glycohemoglobin Standardization Program (NGSP). In addition, a full blood count was also done for all participants, and ethylenediaminetetraacetic acid (EDTA) treated blood samples were stored at −20 degrees Celsius for DNA extraction and analysis.

Genome-Wide DNA Methylation Sequencing
Genomic DNA was extracted from peripheral blood using the Wizard Genomic DNA Purification Kit (Promega, Madison, WI, USA) according to the manufacturer's instructions. At least 2 µg of DNA (concentrations ranging between 70 and 130 ng/µL) with A260/A280 and A260/A230 ratios ≥ 1.8 was shipped frozen on dry ice, as instructed by Arraystar Inc. (Rockville, MD, USA). Methylated DNA immunoprecipitation (MeDIP) was performed by Arraystar Inc. (Rockville, MD, USA) according to Down et al. [39], with minor modifications as follows. About 1 µg of fragmented DNA was prepared for Illumina HiSeq 4000 sequencing as the following steps: (1) end repair of DNA samples with T4 DNA polymerase, Klenow DNA polymerase, and T4 PNK; (2) a single 'A' base was added to the 3' ends with Klenow (exo minus) polymerase; (3) Illumina's genomic adapters were ligated to DNA fragments; (4) DNA fragments were immunoprecipitated by anti-5-methylcytosine antibody (Diagenode); (5) immunoprecipitated DNA fragments were amplified by PCR amplification; (6)~300-600 bp DNA fragments were extracted by gel purification. The completed libraries were quantified by Agilent 2100 Bioanalyzer (Agilent Technologies, Santa Clara, CA, USA). The libraries were denatured with 0.1 M NaOH to generate single-stranded DNA molecules, captured on Illumina flow cell, amplified in situ. The libraries were then sequenced on the Illumina HiSeq 4000 following the TruSeq SBS Kit v5 protocol. The enrichment of DNA immunoprecipitation was analyzed by qPCR using specific methylated sites at H19 locus and non-methylated sites at GAPDH.

MeDIP-Seq Data Analysis
The enrichment of DNA immunoprecipitation was analyzed by qPCR using specific methylated sites at H19 locus and non-methylated sites at GAPDH. Image analysis and base calling were performed using Off-Line Basecaller software (OLB V1.8). After passing a Solexa CHASTITY quality filter, the clean reads were aligned to the human genome (UCSC HG19) using HISAT2 software (V2.1.0). Briefly, individual bases generated from original image files have quality scores, which reflect the probability whether base calling is correct or not. The score is calculated by CHASTITY Formula. The CHASTITY (C) of each base in the short reads is determined by the intensity of four colors (IA, IC, IG, and IT here), and the formula means "the ratio of the highest (IC here) of the four (base type) intensities to the sum of highest two (IC and IG here)." The CHASTITY (C) should be no less than 0.6 in the first 25 bases. Statistically significant MeDIP-enriched regions (peaks) detected by MACS v2 were identified by comparison to input background, using a q-value threshold of 10 −5 . The peaks in samples were annotated by the nearest gene using the newest UCSC RefSeq database. Differentially methylated regions (DMRs) located within gene promoters (TSS − 2000 bp, TSS + 2000 bp) with statistical significance between the two groups were identified by diffReps (Cut-off: log2FC = 1.0, p-value = 10 −4 ).

Gene Ontology (GO) and KEGG Pathway Analysis
The ontology covers three domains, namely biological process, cellular component, and molecular function. Fisher's exact test was used to determine whether there was more overlap between the DE list and the GO annotation list than would be expected by chance. The p value denotes the significance of GO terms enrichment in the DE genes. The lower the p value, the more significant the GO term; a p value ≤ 0.05 was considered significant. Annotation was performed using standard workflow according to http://geneontology.org/. Pathway analysis was done using the Kyoto Encyclopedia of Genes and Genomes (KEGG) database. The p value (EASE score, Fisher's p-value, or hypergeometric p-value) denotes the significance of the pathway correlated to the conditions. The lower the p value is, the more significant the pathway is; a p value ≤ 0.05 was considered significant.