Genome-wide DNA methylation analysis of extreme phenotypes in the identification of novel epigenetic modifications in diabetic retinopathy

Abstract

Background

Aberrant epigenetic modifications such as DNA methylation may contribute to the pathogenesis of DR. We aimed at elucidating the role of novel DNA methylation modifications in diabetic retinopathy (DR) in patients with type 2 diabetes mellitus (T2DM) using an extreme phenotypic design.

Methods/results

Two consecutive studies were conducted. A cross-sectional study using an extreme phenotypic design was conducted to identify rare methylation modifications that might contribute to DR pathogenesis. A 2-year longitudinal nested case–control study was conducted to validate the results and assess whether these novel methylation modifications could be used as biomarkers for predicting DR onset. A large number of differentially methylated CpG sites were identified in the cross-sectional study, and two (cg12869254 and cg04026387) corresponding to known genes were replicated in the longitudinal study. Higher methylation of cg12869254 significantly correlated with macular RNFL thinning in the superior and nasal subregions, and that of cg04026387 correlated with reduced deep capillary plexus VD in the superior and inferior subregions after adjusting for covariates.

Conclusions

Cg12869254 and cg04026387 hypermethylation may complement the known risk factors that contribute to the pathogenesis of DR and as novel biomarkers for disease prediction.

Introduction

Diabetic retinopathy (DR) represents a microvascular end-organ complication of diabetes and is a leading cause of blindness in the working-age population [1]. DR has an insidious onset that delays early detection and treatment [2]. The traditional pathogenesis of DR is largely driven by risk factors such as hyperglycemia, hypertension, and dyslipidemia; however, they explain less than 12% of the risk of developing DR [3]. Growing evidence suggests a clear polygenic risk of developing DR, highlighting the need for new biomarkers for DR. [4]

Aberrant epigenetic modifications such as DNA methylation may contribute to the pathogenesis of DR [5, 6]. Maghbooli et al. [7] noted changes in peripheral blood DNA methylation in DR patients, suggesting the potential value of DNA methylation as a biomarker for DR. However, DNA methylation studies conducted on type 2 diabetes mellitus (T2DM) cohorts based on traditional epigenome-wide association study (EWAS) approaches have yielded unclear results, failing to identify epigenetic modifications reproducibly associated with DR [8,9,10]. The results of these studies were skewed due to the incomplete phenotypes of the patients included. Moreover, the traditional EWAS approaches used in previous studies are inherently limited to common epigenetic modifications, which means that rare changes with large effects will be missed. Extreme phenotypic design is an emerging study design with superior power for detecting rare variants and has been suggested for the screening and identification of DR markers [11,12,13,14,15]. Further, previous studies followed a cross-sectional design and lack of optical coherence tomography (OCT) and OCTA (OCT angiography) biomarkers, which limits causal inferences and interpretations of their results [8,9,10]. Lastly, existing studies did not collect patients’ systemic and ocular factors, specifically renal function and axial length (AL), which are known to be independently associated with DR onset and progression. [16,17,18,19]

In the present study, we hypothesized that the use of extreme phenotypic designs would greatly increase the power to identify candidate sites, thus providing the ability to search for rare DNA methylation modifications. In the discovery stage, a large number of differentially methylated CpG sites (DMSs) in DR were identified from a cross-sectional study. A nested case–control study was conducted to validate the results and longitudinally assess whether aberrant methylation can be used as a biomarker to predict DR. To determine whether these candidate sites play a role in DR pathogenesis, the association between the methylation degree of these sites with macular retinal nerve fiber layer (RNFL) thickness and vessel density (VD) was analyzed. Since diabetic kidney disease often develops simultaneously with DR, we also assessed the potential association of these sites with renal function. [20]

Methods

Study participants

This cross-sectional and longitudinal nested case–control study was conducted at Zhongshan Ophthalmic Center (ZOC) in Guangzhou, China. The study protocol was approved by the institutional review board of ZOC, Sun Yat-sen University (2017KYPJ094). The study procedure followed the principles of the Declaration of Helsinki, and written informed consent was obtained from all the enrolled participants.

Forty-three T2DM patients with DR and 92 T2DM patients without DR at baseline were recruited from 2017 to 2019. All patients were diagnosed with T2DM according to the criteria of the American Diabetes Association [21]. DR was defined according to the revised Early Treatment Diabetic Retinopathy Study (ETDRS) grading scale based on 7-field retinal photography and OCT/A imaging. Patients were excluded from the study if they had (1) severe systemic diseases other than T2DM; (2) cognitive impairment, mental illness, or inability to cooperate with questionnaires or examinations; (3) best-corrected visual acuity (BCVA) less than 20/25; (4) intraocular pressure (IOP) greater than 21 mm Hg; (5) spherical equivalent greater than +3.0 diopters (D) or less than -6.0 D; (6) glaucoma, age-related macular degeneration, amblyopia, and other eye diseases; (7) history of previous intraocular surgery or any other ophthalmic treatment; (8) pterygium and severe refractive media clouding.

Among these patients, 10 were DR-free despite > 20 years of T2DM. All 92 patients without DR at baseline were followed up for 2 years, and 24 of them had new-onset DR. Two consecutive studies were conducted. The cross-sectional study serves as the discovery stage, and the results were then validated in the longitudinal study. The cross-sectional study included (1) a control group (n = 10) diagnosed with T2DM for at least 20 years without any clinical signs or visual-related symptoms of DR, referred to as the NDR group, and (2) a case group (n = 10) of patients diagnosed with DR within 4 years of T2DM diagnosis matched 1:1 by age, sex, and glycated hemoglobin (HbA1c), referred to as the DR group. The longitudinal study enrolled patients without DR at baseline within 2 years of T2DM diagnosis and divided them into two groups according to DR outcome during the 2-year follow-up: (1) DR high-risk group (n = 10): patients with DR onset within 2 years, referred to as the incidence-DR group; and (2) DR stable group (n = 10): patients who remained NDR during the 2-year period matched similarly to the incidence-DR group, referred to as the stable-NDR group (Additional file 2: Fig. S1).

Clinical examinations and SS-OCT/A image acquisition protocol

We obtained a 2 mL fasting venous blood sample anticoagulated with disodium ethylenediaminetetraacetate from each participant and stored at − 80 °C without repeated freezing and thawing. Questionnaires and systemic and ocular examinations were administered to all participants. Details are as follows. After dilatation, seven fundus photographs were obtained according to the ETDRS-7 protocol, using a digital fundus camera (Canon CX-1, Tokyo, Japan). All participants completed a standardized questionnaire to collect information on sex, age, and medical history. Systolic blood pressure, diastolic blood pressure, height (m), and weight (kg) were measured. Body mass index (BMI) was defined as the weight divided by the square of height. Venous blood and urine samples were collected to assay for serum total cholesterol, high-density lipoprotein cholesterol, low-density lipoprotein cholesterol, triglycerides, HbA1c, ACR (albumin creatinine ratio), and estimated glomerular filtration rate (eGFR). The eGFR was calculated using the Cockcroft-Gault formula [22]. The IOP was measured using a non-contact IOP meter (mputeCT-1 Corized Tonometer, Topcon Ltd., Topcon). BCVA was measured using an ETDRS LogMAR visual acuity chart (Precision Vision, Villa Park, Illinois, USA). The AL was measured using a Lenstar LS 900 biometer (HAAG-Streit AG, Köniz, Switzerland).

A commercial SS-OCT imaging system (Triton; Topcon Inc., Tokyo, Japan) was used to obtain structural OCT and OCTA images. The instrument has axial and lateral resolutions of 8 and 20 µm, respectively, and uses a wavelength of 1050 nm with a scan speed of 100,000 A-scans per second. Transverse section images of the fundus and en-face images of the retinal microvascular systems were obtained using a 3D Macula Cube 7 × 7 mm scan mode with a scan density of 512 A-scans × 512 B-scans and Angio Macula 3 × 3 mm scan mode with a scan density of 320 A-scans × 320 B-scans, both centered on the fovea. The transverse section images were automatically segmented into multiple layers, and measurements of the RNFL thickness, including measurements of the superior, inferior, nasal, and temporal regions, were obtained. The en-face images of the retinal vascular systems were automatically segmented into four slabs, including the superficial capillary plexus (SCP) and deep capillary plexus (DCP), and the VDs of each subregion (superior, inferior, nasal, and temporal) in both SCP and DCP slabs were obtained. The results of the automated segmentation were evaluated by an image expert, and the participants’ data were masked during processing. Manual adjustments were performed if segmentation errors were present.

Genome-wide DNA methylation analysis

Whole blood DNA was extracted from fasting venous blood samples using the DNeasy Blood and Tissue Kit (Qiagen, Hilden, Germany). The purity and concentration of the DNA were estimated using a NanoDrop 2000 spectrophotometer (Thermo Scientific). Qualified samples were stored at -20 degrees. Approximately 500 ng of genomic DNA from each sample was used for sodium bisulfite conversion using the EZ DNA methylation Gold Kit (Zymo Research, USA), following the manufacturer’s standard protocol. An Infinium MethylationEPIC BeadChip Kit (Illumina Inc., USA) was used to perform genome-wide DNA methylation analyses. The distribution of cell types was inferred using the robust partial correlation method based on DNA methylation signatures of the constituent cell types.

Probe filtering was performed according to the manufacturer’s protocol to ensure valid analysis. Specifically, all (1) probes with detection p values ≥ 0.01, (2) probes with a bead count less than 3 in more than 5% of samples, (3) non-CpG probes, (4) multi-hit probes, (5) SNP-related probes (within 5 bp of an SNP), and (6) probes located in X and Y chromosomes were filtered out.

The experimental quality of each stage of the genome-wide methylation analysis was carefully assessed using internal reference samples. Staining controls were used to examine the sensitivity and efficiency of the X-staining stage. The efficiency of single-base extension at the X-staining stage was examined using extension controls. Hybridization controls were used to examine the efficiency of DNA hybridization. Target removal controls were used to determine the efficiency of stripping off DNA templates after the single-base extension phase. Bisulfite conversion controls were used to evaluate the conversion efficiency of sulfite to genomic DNA. Specific controls were used to detect non-specific extensions of Infinium I and Infinium II probes. The overall performance of the assay and the sample quality were assessed using non-polymorphic controls. Finally, random sequences that did not hybridize with the target DNA were used as negative controls. The β value was calculated and normalized for quantifying the degree of methylation (Additional file 1: Methods).

Statistical analyses

R software (V4.1.0) was used for all data analyses and the presentation of the results. Continuous variables were presented as means ± standard deviations (SD), and categorical variables were presented as counts (proportions). Continuous variables were compared using Student’s t test, and categorical variables were compared using chi-square tests. To derive methylation levels, array data were analyzed using the ChAMP package (V2.14.0). Batch effects were corrected using the ComBat method. The limma package (V3.40.6) was used to perform linear regression and modified t tests, with age, sex, serum creatinine, AL, and cellular composition included as covariates. Sites that met the significance level with a |Δβ|≥ 0.1 were used for further analyses. Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses were performed to analyze the function of differentiated genes. The Benjamini–Hochberg method was used to reduce the false discovery rate. The association between the degree of methylation of the identified sites and OCT/A measurements, renal function, and metabolic indexes was analyzed using multivariate linear regression. Two-sided p values less than 0.05 were considered statistically significant.

Results

Characteristics of the study groups

The demographics and characteristics of the subjects are summarized in Table 1. In the cross-sectional analysis, subjects from the DR and matched NDR groups shared similar characteristics, with the exception of BCVA (P = 0.007). Similarly, none of the systemic or ocular factors showed any significant differences between the incidence-DR and matched stable-DR groups in the prospective analysis (all P > 0.05). This indicated that the baseline characteristics of the control and case groups were well matched. The baseline characteristics of patients in the incidence-DR group in the prospective analysis were equally comparable to those of patients in the DR group in the cross-sectional analysis (all P > 0.05) (Additional file 2: Table S1).

Table 1 Demographic and ocular characteristics of the study population

Genome-wide DNA methylation analysis

Comparisons revealed significant differences in cellular composition, both in cross-sectional and in longitudinal studies (P < 0.05) (Additional file 2: Figs. S2–S5). The DNA quality was adequate with values of 1.8–2.0 for optical density (OD) 260/OD280. The whole analysis was reliable given that all analytical stages passed rigorous quality control (Additional file 2: Figs. S6–S13). The methylation profiles were compared separately in two comparison groups: DR group versus NDR group, and incidence-DR group versus stable-NDR group. A total of 311 and 413 novel DMSs were identified in the cross-sectional and longitudinal studies, respectively (Fig. 1a, b). The DMSs identified in the cross-sectional study were contained in 141 genes, including CACNA1C, COL19A1, FLT3, S100A13, ZDHHC23, and SLC25A21. The ten most significantly differentiated genes between the NDR and DR groups are listed in Table 2. Compared to the NDR group, 217 sites (69.8%) were hypomethylated and 94 sites (30.2%) were hypermethylated (Fig. 2). Among them, 59 sites were located in the promoter region, 112 in the body region, 7 in the 3′-untranslated region (3′-UTR), 1 in the ExonBnd region, and 32 in the intergenic region (IGR). Concerning CpG positions, 43 sites were located in CpG islands, 183 in the open sea, 10 in the shelf regions, and 75 in the shore regions (Additional file 2: Fig. S14A, B).

Fig. 1

Volcano plot for differential DNA methylation analysis of A DR group and NDR group, and B incidence-DR group and stable-DR group. DR: diabetic retinopathy; NDR: non-DR

Table 2 Top 10 significant genes harboring differentially methylated CpG positions (DMPs) between NDR group and DR group

Fig. 2

Heat map illustrating differentially methylated sites of DR group and NDR group. Red and blue represent levels of hypermethylation and hypomethylation. DR: diabetic retinopathy; NDR: non-DR

The longitudinal nested case–control study identified 197 novel genes, including MICAL1, NLGN1, CACNG2, ZDHHC23, and SLC25A21. The ten most significantly differentiated genes between the stable-NDR and incidence-DR groups are listed in Table 3. Compared to the stable-NDR group, 248 sites (60.0%) were hypomethylated and 165 sites (40.0%) were hypermethylated (Fig. 3). Further analysis of the DNA functional domains revealed that 102 sites were located in the promoter region, 151 in the body region, 6 in the 3’-UTR, 1 in the ExonBnd region, and 153 in the IGR. Regarding CpG positions, 62 sites were located in CpG islands, 219 in the open sea, 15 in the shelf regions, and 117 in the shore regions (Additional file 2: Fig. S14C, D). The results of GO and KEGG pathway analyses are shown in Additional file 2: Figs. S15 and S16.

Table 3 Top 10 significant genes harboring DMP between stable-NDR group and incidence-DR group

Fig. 3

Heat map illustrating differentially methylated sites of incident-DR group and stable-NDR group. Red and blue represent levels of hypermethylation and hypomethylation. DR: diabetic retinopathy; NDR: non-DR

Replicated DMSs in the longitudinal study

The two novel DMSs corresponding to known genes identified in the cross-sectional analysis were replicated in the longitudinal analysis. Cg12869254 and cg04026387, contained in the ZDHHC23 and SLC25A21 genes, were the two novel sites that showed significant methylation changes in the same direction in both cross-sectional (cg12869254: β difference = + 0.11, P = 0.012; cg04026387: β difference = + 0.11, P = 0.018) and longitudinal analyses (cg12869254: β difference = + 0.13, P = 0.015; cg04026387: β difference = + 0.10, P = 0.020). These two sites are both located in the body region concerned with functional genomic distribution and in the open sea relative to CpG islands. Methylation modifications of these two novel DMSs have not been identified using previous traditional GWAS approaches. Hypermethylation at these two sites may serve as a new biomarker for predicting DR.

Association of identified DMSs with OCT/A measurements

In the univariate linear model, the degree of methylation of cg12869254 was negatively correlated with macular RNFL thickness in the superior (β = − 2.6, 95%CI − 5.1, − 0.2; P = 0.039) and nasal subregions (β = − 2.6, 95%CI − 4.4, − 0.7; P = 0.009). The correlation remained stable in the multivariate model after adjusting for age, sex, T2DM duration, and HbA1c level (β = − 3.0, 95%CI − 5.8, − 0.2; P = 0.035 and β = − 2.2, 95%CI − 4.4, − 0.1; P = 0.049) (Table 4). The degree of methylation of cg04026387 was found to be negatively correlated with macular DCP VD in the superior (β = − 0.8, 95% CI − 1.6, − 0.1; P = 0.049) and inferior subregions (β = − 1.4, 95% CI − 2.4, − 0.4; P = 0.008) and positively correlated with SCP VD in the nasal subregion (β = 1.9, 95% CI 0.5, 3.3; P = 0.011) in the univariate linear model. The correlation remained stable in model adjusting for age, sex, T2DM duration, and HbA1c (β = − 0.8, 95% CI − 1.7, − 0.1; P = 0.039, β = − 1.5, 95% CI − 2.8, − 0.2; P = 0.028, and β = 2.2, 95% CI 1.1, 3.3; P = 0.001) (Table 5).

Table 4 Association between cg12869254 with OCT/A measurements in univariable and multivariable linear regression models

Table 5 Association between cg04026387 with OCT/A measurements in univariable and multivariable linear regression models

Association of identified DMSs with renal function and metabolic indexes

The degree of methylation of cg12869254 was significantly correlated with reduced eGFR (β = − 0.2, 95% CI − 0.4,-0.01; P = 0.038), and that of cg04026387 was significantly correlated with elevated ACR (β = 0.7, 95% CI 0.2, 1.1; P = 0.014), after adjusting for age, sex, T2DM duration, and HbA1c level (Additional file 2: Tables S2 and S3). Several metabolic indexes were also related to the methylation of these two sites (Additional file 2: Table S4, Additional file 1: Results).

Discussion

This study consisted of genome-wide DNA methylation analyses of extreme phenotypes within two consecutive studies that identified novel epigenetic modifications of DR. In the discovery stage, a large number of novel candidate sites were identified. To validate the results and assess whether these methylation modifications could be used as biomarkers to predict DR, we analyzed DNA methylation in an independent longitudinal study. Our results suggest that aberrant methylation may predict DR, as hypermethylation at cg04026387 and cg12869254 identified in the cross-sectional study was successfully replicated in the longitudinal study. Finally, the methylation degree of these two sites negatively correlated with macular RNFL thickness and DCP VD, which further confirmed their potential for indicating DR [23, 24]. Therefore, methylation at these two sites may contribute to neural and vascular damage in the early stages of DR and serve as novel biomarkers for predicting the disease.

In this study, an extreme phenotypic design was adopted to identify the specific DMSs, which enabled us to identify novel methylation modifications, despite the relatively small yet sufficient sample size, as in studies of similar design [25,26,27,28]. Conversely, the results of previous studies on DR were skewed by the incomplete phenotypes of included participants, that is, healthy individuals and/or inclusion of individuals with diabetes diagnosis of < 20 years as DR control groups, and/or DR patients with random DR-free T2DM duration as the DR case group. Current evidence suggests that an extreme phenotypic design is necessary for screening and identifying DR markers [29, 30]. The extreme phenotypic design is a recent concept that is used to increase homogeneity and better distinguish signals from noise in genome-wide analysis. Under this approach, the decreased variation in homogeneous samples results in a robust ability to identify signals, even in a relatively small number of samples. Concurrently, it enables the discovery of rare methylation alternations, which is another advantage of this design. Thus, we included patients diagnosed with T2DM for at least 20 years without any signs of DR as an extreme control group, and patients diagnosed with DR within 4 years of T2DM diagnosis were included as an extreme case group. As most patients typically develop DR within 17 years of T2DM diagnosis [11], while epigenetic changes typically initiate retinal impairment within a short period, the enrolled patients may represent extreme phenotypes enriched in “protective” and “risk” epigenetic alternations, respectively. Therefore, this model increased the candidate site identification in this study.

To our knowledge, this is the first study to longitudinally investigate the relationship between DNA methylation and T2DM in patients with DR. Thus far, all previous studies on T2DM patients have been cross-sectional, which limits causal inferences of their results. However, DNA methylation has been suggested to be dynamic rather than rigidly fixed in various diseases and aging process [31]. Therefore, it is important to clarify the time-ordered relationships between such aberrant methylation and DR pathogenesis. Our cross-sectional study revealed significant hypomethylation in sites within S100A13 in the gene-positive DR group, which is in agreement with the results obtained by Li et al. [8] Demethylation of this gene purportedly increases hyperglycemia-induced damage through calcium signaling and the RAGE pathway. However, these results were not replicated in the longitudinal study. Similarly, many other candidate sites identified in the cross-sectional discovery stage exhibited no significant differences in the longitudinal study. This highlights the limitations of previous cross-sectional studies in which the chronological relationship between exposure and outcome was unclear. The changes in methylation at these sites may be explained as concomitant consequences of DR development rather than as primary contributing factors to DR pathogenesis or progression. Cg12869254 and cg04026387 were successfully replicated in this longitudinal study, suggesting that methylation alterations at these two sites precede DR, supporting their role in the early pathogenesis of DR.

The two novel sites identified in this study are included in the ZDHHC23 and SLC25A21 genes, which have not been identified in previous studies. A recent animal study demonstrated the neurotoxic effects mediated by the downregulation of the ZDHHC23 gene [32]. Decreased ZDHHC23 levels dysregulated interactions between palmitoyl-protein thioesterase-1 and acyl protein thioesterase-1, leading to increased plasma membrane H-Ras and subsequent microglial proliferation in mouse brains. Activation of microglia leads to increased secretion of tumor necrosis factor-alpha, interleukin-6, monocyte chemoattractant protein-1, and complement component C1q, promoting the transformation of astrocytes to the neurotoxic A1 phenotype [33]. The astrocytes of the A1 phenotype not only lose their original neurotrophic role, but also actively produce neurotoxins that degrade neurons. Retinal neurodegeneration is considered an early component of DR, which can precede visible vasculopathy [23, 24]. In this study, the degree of methylation of cg12869254 negatively correlated with macular RNFL thickness in the superior and nasal subregions, while the RNFL thickness represented the axons of retinal ganglion cells. Our teams and other groups have demonstrated that reduced RNFL thickness was significantly associated with a higher risk of development and progression of DR [34]. Therefore, we hypothesize that ZDHHC23 gene-mediated neurotoxicity caused by hypermethylation of cg12869254 may also contribute to the retinal neurodegeneration in DR. Additional studies are warranted to further elucidate the potential role of ZDHHC23 in DR.

SLC25A21 is involved in intracellular organic acid catabolism and metabolism [35]. Within mitochondria, it transports 2-oxoadipate and 2-oxoglutarate to be converted into acetyl-CoA. Boczonadi et al. [36] found that the dysfunction of SLC25A21 may lead to impaired transport of 2-oxoadipate, which in turn disrupts the degradation of lysine and tryptophan, leading to the accumulation of 2-oxoadipate, piperidinic acid, and quinolinic acid. Overexposure to these cytotoxic substances leads to irreversible mitochondrial and tissue damage mediated by free radicals. Consistently, the methylation degree of cg04026387 locus negatively correlated with macular DCP VD in our study, suggesting that it is a more sensitive indicator than VD in other plexus [37]. Therefore, we hypothesize that similar mechanisms of mitochondrial damage and vascular damage may also occur in the pathogenesis of early DR. However, SCP VD in the nasal subregion appears to be increased with higher methylation of cg04026387. Similar results have been reciprocated in studies on rats and humans [38, 39]. One possible explanation for this is that the deeper capillary layer is more sensitive to mitochondrial damage. In the early stages of DR, the above-mentioned factors may act mainly on the DCP, whereas the SCP exhibits a compensatory flow increase secondary to the reduced DCP VD. However, further studies investigating these genes are required to validate our hypothesis.

DKD is another major microvascular complication of diabetes mellitus [40]. Given the commonalities in pathogenesis and frequent co-development of DR and DKD, we further analyzed the relationship between these two sites and renal function indexes [20, 41,42,43]. Hypermethylation at these two identified sites was significantly associated with elevated ACR and reduced eGFR, further increasing our confidence in the involvement of these two sites in diabetic microvascular complications. We hypothesize that inflammatory responses, mitochondrial dysfunction, and oxidative stress discussed above may also damage the endothelial cells of glomerular capillaries [41, 43,44,45], thus leading to urinary protein leakage and decreased glomerular filtration function.

To our knowledge, few studies have collected systemic and ocular factors, such as renal function and AL, while analyzing the molecular etiology of DR. These factors have been independently associated with DR onset and progression [16,17,18,19]. Thus, the differentiated sites identified in previous studies might have been confounded by the presence of these factors. This study obtained various systemic and ocular factors from the participants, and most characteristics were well matched in each comparison group. This increased the reliability of our results.

This study had certain limitations. First, since most patients typically develop DR within 17 years of T2DM diagnosis [11], patients diagnosed with T2DM for 20 years without DR were included in the gene-negative control (NDR) group. However, even if the odds are low, these patients may still develop DR in the future. Future longitudinal studies are needed to determine whether these protective factors are still effective when the diabetes duration is over 20 years. Second, the sample size of this study was relatively small, and therefore the results should be treated with care. However, the adopted extreme phenotypic design greatly increases genetic homogeneity, which provides a robust ability to distinguish signals from noise. Despite the relatively small sample size, two sites were successfully replicated in the longitudinal study and were correlated with OCT/A parameters.

Conclusion

In conclusion, using extreme phenotypes, this study identified a large number of novel methylation modifications of DR in the discovery stage from a cross-sectional study and further validated these results via a longitudinal study. Cg12869254 and cg04026387 were the two sites corresponding to known genes that were successfully replicated, and the degree of methylation at these sites was associated with macular RNFL and VD. These two novel sites may serve as complementary to the known risk factors that contribute to the pathogenesis of DR and as biomarkers to identify diabetic patients at high risk of incident DR. Monitoring these sites may facilitate risk stratification of the diabetic population and offer a potential tool for screening and therapeutic interventions.