Genome-wide analysis of DNA methylation identifies S100A13 as an epigenetic biomarker in individuals with chronic (≥ 30 years) type 2 diabetes without diabetic retinopathy

Background This study aimed to determine the epigenetic biomarkers of diabetic retinopathy (DR) in subjects with type 2 diabetes mellitus (T2DM). This retrospective study is based on the Shanghai Xinjing community prevention and treatment administrative system of chronic diseases. The subjects enrolled herein were T2DM patients who had undergone long-term follow-up evaluation in the system. Two consecutive studies were conducted. In the discovery cohort, among 19 subjects who had developed DR with a DM duration < 3 years and 21 subjects without DR > 30 years after being diagnosed with DM, an Infinium Human Methylation 850 Beadchip was used to identify differential methylation regions (DMRs) and differential methylation sites (DMSs). The function of the genes was assessed through KEGG enrichment analysis, Gene Ontology (GO) analysis, and pathway network analysis. In the replication cohort, 87 DR patients with a short DM duration and 89 patients without DR over a DM duration > 20 years were compared to assess the association between DMSs and DR upon pyrosequencing. Results A total of 34 DMRs were identified. Genes containing DMSs with the top 5 highest beta value differences between DR and non-DR participants were located on chromosome 1 and were present in the S100A13 gene, which was associated with 71 GO terms. Two S100A13 gene sites, i.e., cg02873163 and cg11343894, displayed a good correlation with DR on pyrosequencing. Conclusions DMSs in the S100A13 gene may be potential biomarkers of DR.


Background
Epigenetic determinants reportedly contribute to the etiology of numerous systemic disorders such as diabetes mellitus (DM) [1,2]. Furthermore, hyperglycemia can affect DNA methylase activity or the methylation level of specific genes, which is the primary epigenetic modification [3][4][5][6][7]. Diabetic retinopathy (DR) is a common complication of DM and a leading cause of blindness among individuals aged 20-62 years in developed countries [8]. However, the exact etiology of DR remains unclear. Chen et al. reported that approximately 300 differential methylation sites (DMSs) in DNA were associated with DR onset in type 1 DM (T1DM) patients [9]. Agardh et al. reported an association between DNA methylation and proliferative DR among individuals with T1DM [5]. However, among individuals with T2DM, no specific gene methylation levels have been reported to be associated with DR, even in PubMed.
In the present study, we assessed the epigenetic biomarkers of DR in T2DM patients. In the discovery cohort, among 19 subjects who developed DR with a DM duration of < 3 years and 21 subjects without DR > 30 years after being diagnosed with DM, > 850,000 gene sites were assessed with a methylation chip. The functions of genes harboring DMSs with different degrees of methylation between the DR and non-DR participants were assessed through KEGG enrichment analysis, Gene Ontology (GO) analysis, and pathway network analysis and assessed for their association with DR. In the replication cohort, 87 DR patients with a short DM duration and 89 patients without DR over a DM duration > 20 years were compared to validate the association between the DMSs identified in the discovery cohort and DR through pyrosequencing.

Discovery study
The basic and clinical characteristics of the 19 DR and 21 non-DR subjects are summarized in Table 1. Except for differences in the duration and hemoglobin A1C (HbA1c) levels, the two groups did not significantly differ in age, sex, glucose levels, body mass index (BMI), spherical equivalent rate, intraocular pressure, or axial length. In the DR group, one had proliferative DR, 5 patients had severe nonproliferative DR, and 13 patients had moderate nonproliferative DR. Diabetic macular edema was detected in 3 patients. The 21 non-DR subjects did not present signs of DR on annual eye examination after they were diagnosed with T2DM. The duration of DM among these subjects was > 30 years.
We identified thirty-four differential methylation regions (DMRs) when comparing the DR group and non-DR group. The DMRs comprised 11 hypermethylated regions and 23 hypomethylated regions. Of the 12 generelated sites, two DMSs were located at the transcriptional start site (TSS) 1500, one at TSS 200, one in the 5′-untranslated region UTR-5, three in EXON1, two in the gene body, and three in the 3′-UTR (Table 2).
Global DNA methylation levels did not significantly differ between the DR group and the non-DR group (P > 0.05). Furthermore, among 290 DMSs, methylation levels were decreased at 188 sites and increased at the other 102 sites in the DR group (Supplement 1). These DMSs were evenly distributed on all chromosomes (Fig.  1). The genes containing DMSs with the top 5 differences in beta values between DR and non-DR participants are listed in Table 3; these DMSs were detected in the S100A13 gene.
On GO enrichment analysis and KEGG pathway analysis (Supplement 2 and 3, respectively), the top 30 items are summarized in Fig. 2 a and b, respectively. DR-associated pathways included the Wnt signaling pathway, ubiquitin-mediated proteolysis, Tolllike receptor signaling pathway, tight junctions, the transforming growth factor-beta signaling pathway, the phosphatidylinositol signaling pathway, natural killer cell-mediated cytotoxicity, inositol phosphorylation metabolism, the hedgehog signaling pathway, aminoglycan degradation, Fcγ receptor-mediated autophagy, cell adhesion factor signaling, the calcium signaling pathway, regulation of muscle contraction, neural tube formation, neural cell development, iron ion balance, fibroblast proliferation, phagocytosis associated with Fcγ receptor pathogenesis, glucosestimulated cellular responses, and developmental processes. Differential methylation in apoptotic responses and genes involved in AMPK activity are indicated in Fig. 2 a, b, and c. S100A13 was associated with 71 GO terms, including signal transduction and calcium ion binding, which are reportedly associated with DR pathogenesis [10,11]. Therefore, the DMSs of the S100A13 gene were considered candidate DMSs.

Replication study
The basic and clinical characteristics of the 87 DR and 89 non-DR subjects are summarized in Table 4. Except for DM duration and HbA1C levels, no significant differences were observed in blood glucose levels, BMI, intraocular pressure, equivalent spherical rate, and axial length between the two groups. The 89 participants in the non-DR group did not present signs of DR on annual eye examination after being diagnosed with T2DM. Their disease durations were > 20 years. The association between DR and DMS in the S100A13 gene was verified through pyrosequencing (primer sequences are provided in Supplement 4). Significant differences (P < 0.05) between the DR and non-DR groups were observed at sites cg02873163+12 (41.91 ± 12.60%; 47.73 ± 12.36%), cg02873163+2 (54.99 ± 14.77%; 61.45 ± 13.68%), cg02873163 (51.03 ± 14.63%; 57.21 ± 13.26%), cg11343894+5 (45.86 ± 15.88%; 51.47 ± 14.34), and cg11343894 (48.83 ± 15.17%; 54.51 ± 14.45) in the S100A13 gene. After adjusting for age, sex, BMI, HbA1C, and mean arterial pressure, the average methylation level of all 5 sites of cg02873163 and cg11343894 was positively associated with DR through binomial logistic regression analysis (Table 5).

Discussion
To identify possible epigenetic biomarkers of DR among individuals with T2DM, we enrolled T2DM patients segregated into a discovery cohort and a replication cohort and identified 34 DMRs associated with DR in Chinese individuals with T2DM. These 34 DMRs included 23 (67.65%) hypomethylation sites, concurrent with a previous report on T1DM patients [5]. DR results from a cascade of enzymatic reactions. With the regulation of DNA methyltransferases and ten-eleven-translocation proteins, DNA hypermethylation/demethylation alters the binding of transcription factors and genes in most cases [12][13][14]. Interactions between genetic and epigenetic factors may play a role in DR pathogenesis. Fig. 1 The distribution of differential methylation sites (DMSs). Red markers represent the hypermethylation sites of the DR group, and green markers represent the hypomethylation site Table 3 Genes harboring differential methylation sites (DMSs) with the top 5 differences in β values between diabetic retinopathy (DR) and non-DR subjects in the discovery study  Chen et al. reported that different blood glucose levels affect the extent of DNA methylation in specific genes in whole blood or monocytes, which subsequently affects DR pathogenesis [9]. Herein, blood glucose and HbA1c levels were lower among DR subjects than among non-DR subjects, indicating that the etiology of DR is complex and does not simply depend on glucose levels. We show that the S100A13 gene is potentially involved in DR pathogenesis and that the cg02873163 and cg11343894 sites in the S100A13 gene may serve as potential biomarkers of DR.
The S100A13 gene is involved in calcium ion homeostasis, energy metabolism, inflammation, apoptosis, regulation of proliferation, differentiation, and interactions with transcription factors and nucleic acids within cells and activates surface receptors, including the receptor for advanced glycation end-products and Toll-like receptor 4, G-protein-coupled receptors, scavenger receptors, or heparan sulfate proteoglycans and N-glycans [15,16]. Furthermore, S100A13 regulates calcium levels and fibroblast growth factor-1 and interleukin-1 alpha secretion through a noncanonical pathway [10,11,17]. Mandinov et al. reported that S100A13 mediates IL-6induced damage to the macrovasculature [18]. We speculate that significant demethylation of S100A13 in individuals with DR downregulates S100A13 and upregulates p38 MAPK and nuclear factor-kappa B through calcium signaling and the RAGE pathway, thus increasing the damage caused by hyperglycemia [17,[19][20][21][22][23]. Patients harboring hypomethylated S100A13 may require closer follow-up evaluation and more stringent blood glucose screening.
It is generally difficult to confirm a gene-negative contrast group in genetic studies on chronic diseases. For example, DR may occur throughout the lifetime of an individual with T2DM, and the exact diagnosis of non-DR can only be made towards the end of an individual's lifespan. Herein, we selected non-DR individuals with a maximal DM duration to serve as non-DR participants. Furthermore, we selected individuals who developed DR with a DM duration of < 3 years as the gene-positive group, since we believe that epigenetic changes initiate retinal neural and vascular impairment within a short period. This study cohort displayed large differences in age, onset age, and duration of DM between the two groups, subsequently leading us to speculate whether differences in DNA methylation are associated with age or duration of T2DM. Therefore, we conducted two consecutive studies, including different non-DR contrast groups with different DM durations, and found that (See figure on previous page.) Fig. 2 The GO enrichment map of differential methylation sites (DMSs), and the top 30 items are displayed in the graph. The larger the enrichment value is, the more significant the enrichment is. Fig. 2 b shows the KEGG pathway analysis of DMSs in genes, and the top 30 items are displayed in the graph. Figure 2 c shows the potentially most important pathway network, and the differentially methylated genes are involved in pathways of cytokine receptor signaling, MAPK signaling pathway, apoptosis, calcium signaling pathway, adhesion factor, ErbB signaling pathway, cell cycle, Wnt signaling pathway, mTOR pathway, and others. Figure 2    DMSs in the S100A13 gene were present in both studies; moreover, the associations between DMSs and DR were confirmed with regression analysis. This study has several limitations. First, subjects developing DR within a relatively short duration or those without DR in the long term may harbor gene mutations, which potentially affect approximately 14% of the methylation levels at corresponding or adjacent sites [24]. Second, each cell type has a different DNA methylation programming system [18]. DMSs in whole blood cells can only reflect the potential for the occurrence of diabetic vascular disease instead of duplicate DMSs in retinal cells. Finally, 22 DMRs are present in CpG islands located distal to any known gene among the 34 DMRs, and their functions warrant further elucidation in future studies.

Conclusions
Our study is the first to detect DMRs and DMSs associated with DR in T2DM patients, indicating that DMSs in the S100A13 gene serve as potential biomarkers of DR. Future studies with more genetically diverse cohorts having a shorter DM duration of < 10 years and prospective cohort studies are required to confirm the effect of epigenetic changes in the S100A13 gene on DR incidence.

Methods
This retrospective study is based on the Shanghai Xinjing community prevention and treatment system of chronic diseases, which was established in 1999. The system has curated annually updated health data on all local residents diagnosed with T2DM. In 2003, our group participated in this system, and we have been conducting and supervising studies on the prevention and treatment of ocular diseases [25][26][27][28][29]. From June to August 2016, we conducted an epidemiological study on eye complications of DM upon receiving ethical approval by the Ethics Committee of Shanghai General Hospital, Shanghai Jiaotong University School of Medicine (2013KY023). This study describes the epigenetic analysis conducted within that study. This study adheres to the tenets of the Helsinki Declaration, and the study subjects provided written informed consent to participate in the study.
The inclusion criteria were as follows: (1) Chinese Han background, (2) a diagnosis of T2DM in accordance with the World Health Organization criteria [30], and (3) the ability to comply with all the required examinations. The exclusion criteria were as follows: (1) occurrence of eyelid diseases, strabismus, corneal diseases, lens diseases, and other eye diseases potentially affecting the outcomes of retinal examinations; (2) occurrence of other eye diseases, including glaucoma or macular degeneration, which may cause fundus retinal microvasculopathy; (3) history of eye surgery or trauma; (4) occurrence of severe systemic diseases, including those involving the respiratory system, circulatory system, and excretory system; and (5) history of cancer.
The research staff were fully trained and experienced with epidemiologic studies [25][26][27]. We gathered data including demographic characteristics and conducted a routine eye examination. DR was diagnosed on the basis of the well-accepted international diagnostic criteria. For each participant, 2 mL samples of fasting blood were obtained, the cell composition was determined for further data analysis, and DNA was extracted using a QIAN amp blood kit (Qiagen, Hilden, Germany) in accordance with the manufacturer's instructions.

Discovery study Enrolled participants
The discovery cohort included 19 subjects who developed DR with a DM duration of < 3 years and 21 subjects without DR > 30 years after being diagnosed with DM.

DNA methylation quantification
The genomic DNA (500 ng) was subjected to bisulfite conversion using an EZ DNA methylation kit (Zymo Research, Orange, CA, USA) in accordance with the manufacturer's instructions. An Infinium Human Methylation 850 Beadchip methylation chip (Illumina, San Diego, CA, USA), which covers 99% of all RefSeq genes and contains 867,531 sites, was used. Methylation analysis involved data quality control, preprocessing, methylation difference site analysis, and methylation difference region analysis.
The quality control step for the analysis of the 850k chip included checking for sex (X,Y-chromosome data were not included in the DMS/DMR analysis) [31], filtering the SNP sites (minor allele frequency > 5% and probes with SNPs of minor allele frequency > 5% within 10 base pairs of the CpG sites were excluded) [32], (list of CpG sites is available at http://genetics. emory.edu/research/conneely/annotation-r-code.html), normalization (Subset-quantile Within Array Normalization was used to normalize the DNA methylation probe intensity) [33], following the criteria for retaining probes (detection P < 0.05 in more than half of the samples) and individuals (detection P value of more than 95% of the probes of each sample < 0.05).

Replication study Enrolled participants
Pyrosequencing was performed to verify the correlation between the candidate DMSs and DR. The replication cohort included 87 subjects who developed DR with a DM duration of < 3 years and 89 subjects without DR over a DM duration > 20 years.
Pyrosequencing DNA (500 ng) was subjected to bisulfite conversion with the EZ DNA methylation kit (Zymo Research) in accordance with the manufacturer's instructions. Pyrosequencing was performed with the PyroMark Q96 ID system (Qiagen, Valencia, CA, USA) in accordance with the manufacturer's recommendations. The PyroMark PCR Master Mix kit (Qiagen), streptavidin-coated beads (GE Healthcare, Uppsala, Sweden), PyroMark Gold Q96 reagents (Qiagen), PyroMark Q96 Vacuum Workstation, and PyroMark Q96 software (version 2.5.8, Qiagen) were used to determine and analyze DMSs.

Data processing and statistical analysis
The clinical data were statistically analyzed using the SPSS 22.0 software (IBM Corporation, Armonk, NY, USA). The data are expressed as the mean ± standard deviation values, and calculated values are presented as frequencies and percentages. Data normality was assessed using the Kolmogorov-Smirnov test. Normally distributed data were analyzed using Student's t-test to compare differences in various values between the DR and non-DR groups. Skewed data were analyzed using the Mann-Whitney test. Statistical tests performed herein were the chi-square test and independent samples t test.
Minfi software package (version 1.25.1, http://www. bioconductor.org/packages/release/bioc/html/minfi. html) and IMA software package (version 2.0, http:// www.rforge.net/IMA), both written with R software, were used to determine the differences in the degree of DNA methylation between the DR participants and the non-DR participants [34]. Minfi package was used to calculate the beta value (interval between 0 and 1), an indicator of the degree of methylation using the following equation: beta i ¼ maxðy ði;methyÞ ;0Þ maxðy ði;methyÞ ;0Þþ maxðy ði;unmethyÞ ;0Þþ100 . The distribution of the beta values was skewed, and all outlier beta values were excluded (omitted values > 4 SD).
The sitetest algorithm and the regionswrapper algorithm, both from the IMA package, were used to define DMS and DMR, respectively. Du et al. suggested the threshold of beta value as between 0.05 and 0.15 [35], and according to the Illumina's instructions, 0.14 is the minimum value recommended to guarantee the proper sensitivity. Therefore, if differences between beta values at a certain region or site between the DR and non-DR groups were greater than 0.14 and the P value determined via a pooled t test was less than 0.05, the region was considered a DMR, and the site was considered a DMS without regression analysis.
The functions of genes containing DMRs or DMSs were assessed via KEGG enrichment analysis, GO analysis, and pathway network analysis. We mapped the methylated sites to GO terms in the Gene Ontology database and determined and corrected (Bonferroni) the number of genes for each GO term. KEGG pathway analysis (http://www.kegg.jp/) was performed to confirm the enriched pathways. The threshold for statistical significance (< 0.05) was defined on the basis of the corrected P value. The gene coexpression network was constructed in accordance with the normalized signal intensity to identify gene interactions. A strong correlation (Pearson correlation > 0.9) indicated the interaction between two genes. Cytoscape (National Institute of General Medical Sciences, Boston, MA, USA) was used to construct coexpression networks. If a gene harbored a DMS with higher differences in beta values between DR and non-DR participants and the corresponding protein was involved in previously reported pathways associated with DR pathogenesis, the DMS was considered a candidate DMS in the following replication study. The R 3.5.1 pheatmap package was used to calculate and draw the heat map.
In the replication part, the association between methylation level (%) on pyrosequencing and DR was adjusted for demographic characteristics, including age, sex, Hb1Ac (%), body mass index, and mean arterial pressure, which were assessed through binomial logistic regression analysis. A P value < 0.05 was considered statistically significant.