Differential Methylation of Telomere-Related Genes Is Associated with Kidney Disease in Individuals with Type 1 Diabetes

Diabetic kidney disease (DKD) represents a major global health problem. Accelerated ageing is a key feature of DKD and, therefore, characteristics of accelerated ageing may provide useful biomarkers or therapeutic targets. Harnessing multi-omics, features affecting telomere biology and any associated methylome dysregulation in DKD were explored. Genotype data for nuclear genome polymorphisms in telomere-related genes were extracted from genome-wide case–control association data (n = 823 DKD/903 controls; n = 247 end-stage kidney disease (ESKD)/1479 controls). Telomere length was established using quantitative polymerase chain reaction. Quantitative methylation values for 1091 CpG sites in telomere-related genes were extracted from epigenome-wide case–control association data (n = 150 DKD/100 controls). Telomere length was significantly shorter in older age groups (p = 7.6 × 10−6). Telomere length was also significantly reduced (p = 6.6 × 10−5) in DKD versus control individuals, with significance remaining after covariate adjustment (p = 0.028). DKD and ESKD were nominally associated with telomere-related genetic variation, with Mendelian randomisation highlighting no significant association between genetically predicted telomere length and kidney disease. A total of 496 CpG sites in 212 genes reached epigenome-wide significance (p ≤ 10−8) for DKD association, and 412 CpG sites in 193 genes for ESKD. Functional prediction revealed differentially methylated genes were enriched for Wnt signalling involvement. Harnessing previously published RNA-sequencing datasets, potential targets where epigenetic dysregulation may result in altered gene expression were revealed, useful as potential diagnostic and therapeutic targets for intervention.


Introduction
The incidence of diabetes is increasing worldwide, with cases expected to rise to 578 million by 2030 [1]. Approximately 30 to 40% of individuals with diabetes go on to develop kidney disease [2]. In 2018, diabetes was the most common primary disease in patients requiring renal replacement therapy (RRT) in the UK and over 40% of incident endstage kidney disease (ESKD) patients in the USA had diabetes, with significant associated healthcare costs [3,4].

Materials and Method
A flowchart summarising the experimental design is shown in Figure 1.

Study Populations
All participants provided written informed consent and were recruited from one of four populations defined in Supplementary materials.

Materials and Method
A flowchart summarising the experimental design is shown in Figure 1.

Study Populations
All participants provided written informed consent and were recruited from one of four populations defined in Supplementary materials.

Telomere Length
Telomere length values were determined via monochrome quantitative polymerase chain reaction (qPCR) using relative ratio of telomere repeat copy number to a single copy gene 36B4 [60] in DNA derived from whole blood. Further details are included in Supplementary materials.

Mendelian Randomisation
Two-sample Mendelian randomisation was performed as described in Supplementary materials, harnessing SNPs utilised by Park et al. [50] and Codd et al. [44].

Epigenome-Wide Association Study (EWAS)
Existing quantitative (blood-derived) DNA methylation levels were extracted for genes relevant to telomere function from the UK-ROI cohort [63,64]. Briefly, 485,577 methylation sites for 150 DKD cases were compared to 100 controls with T1D and no evidence of kidney disease; data from 1091 CpG sites in 376 telomere-related genes were extracted for this analysis. Venn diagrams were created using the ggVenndiagram package (v. 1.2.0) in R Studio (v. 1.4.1717) [65,66].

Fine Mapping of Top Methylation Sites
Following bisulphite conversion using the EZ DNA Methylation-Lightning™ (Zymo Research, Freiburg im Breisgau, Germany), fine mapping was performed for 7 of the top-ranked differential methylation sites in MAD1L1. Primer pair sequences and PCR conditions are shown in Table S2. Thermal cycling was performed on the MJ Tetrad PTC-225 thermal cycler with the following conditions: 15 min at 95 • C; 35 cycles of 45 s at 94 • C, 45 s at 55-65 • C (temperature optimized for each primer pair) and 1 min at 72 • C; 10 min at 72 • C. Samples were subjected to exonuclease and phosphatase clean-up and sequencing reaction using BigDye Terminator v3.1 Ready Reaction Cycle Sequencing kit (Applied Biosystems, Massachusetts, MA, USA) before cycle sequencing on the ABI 3730 ® capillary sequencer. Results were evaluated using Contig Express on Vector NTI version 11.5.1 (Life Technologies, California, CA, USA).

Gene Ontology Analysis
Gene ontology analysis of biological processes, with subsequent clustering of gene ontology terms, was performed using the ViSEAGO package (v. 1.6.0), in R studio (v. 1.4.1717) [66,67], according to processes described in Supplementary materials.  [69] and the GWAS catalogue [70] were also harnessed, searching within 150 kilobases either side of the SNP co-ordinates.

Differential Expression Analysis
Gene expression data from previous RNA-sequencing analyses of DKD and control tissues [71,72] were analysed as described in Supplementary materials.

Statistical Analysis
Plots in Figure 2 were created using R packages ggplot2 (v. 3.3.5) and ggpubr (v 0.4.0) [73,74]. These boxplots display the median, together with the first and third quartiles (hinges), with whiskers extending from the hinges to values at most ±1.5* interquartile range. When comparing relative telomere length (RTL) between groups, Kruskal-Wallis and Wilcoxon rank sum tests were performed (using ggpubr compare_means functions) due to the non-normality of RTL data. Where relevant, p-values for these tests were adjusted using the Bonferroni correction method. Where corrected RTL values are shown, RTL values were corrected for the covariates age, duration of T1D, sex, body mass index (BMI) and glycated haemoglobin (HbA1c) via a series of linear regressions, harnessing residuals as the corrected RTL values.  Figure 2B) (n = 536 cases/611 controls, p = 6.6 × 10 −5 ). This association was ameliorated when RTL was corrected for age, T1D duration, sex, BMI and HbA1c (p = 0.028). When stratifying by age group, a significant difference in RTL in individuals with DKD, compared to control individuals, was only observed in the 46 to 55 age group (p = 0.0087) (Figure 2A). In the Steno replication cohort, however, no significant difference in telomere length (assessed by qPCR) was observed in DKD compared to controls (n = 78 cases/153 controls, p = 0.3). RTL in 157 cases and 116 controls from this population were previously measured using Southern blot [56]. The correlation between the qPCR and Southern blot (SB) results (r > 0.8) highlighted the consistency of methods used [75,76].

Genetic Variants Were Nominally Associated with Telomere Length
Existing genome-wide SNP data were used to investigate association with RTL in 1019 individuals from the UK-ROI cohort, adjusting for age ( Figure S1). Nominally associated hits (p < 0.005) are shown in Table S3. Of the SNPs within the 376 telomere-related genes (Table S1), no SNP reached significance (p = 10 −8 ); however, 36 SNPs in 22 genes demonstrated nominal association (Table 1 and Table S4). In total, 17 of these SNPs were within genes nominally associated with hypertension, diabetes, renal, glomerulonephritis or nephritis-related phenotypes in the UK Biobank (Table S4). Moreover, searching within 150 kb of these SNP co-ordinates revealed that 25 of these regions (69.4%) were significantly associated with a range of renal or cardiovascular phenotypes in the Common Metabolic Disease portal database (clumped by linkage disequilibrium) [69], which aggregates and analyses human genetic and functional genomic information linked to common metabolic diseases from up to 398 datasets (Table S5). Using these same regions to search the GWAS Catalogue [70], it was shown that the rs852540 SNP region (7:5,383,963-5,733,963) contained the variant rs7808152, previously associated with telomere length (p = 1 × 10 −6 , Beta = −0.1602507, CI = 0.096-0.225) in a cohort of 902 European-ancestry individuals (Netherlands) [77]. It was also shown that, within 150 kb of variants rs2209437, rs2025557 and rs1536078 (all with the closest gene SH3GL2), variants associated with DNA methylation (rs7032102, p = 2 × 10 −8 ) and epigenetic clock age acceleration (GrimAge) (rs1114790, p = 10 −8 , Beta = −1.0232, CI = 0.67-1.38) were identified in a cohort of up to 954 individuals from the UK [78]. The highest-ranked telomere-related SNP was rs2725385 in the WRN gene (p = 2.09 × 10 −4 , OR = 1.28, 95% CI = 1. 22-1.45). An additional 5 of the top 36 SNPs were also located in this gene, all within strong linkage disequilibrium (D > 0.8).

Genetic Variants within Telomere-Related Genes Were Nominally Associated with DKD and ESKD
Using existing GWAS data, SNPs within the telomere-related genes were investigated for an association with DKD and ESKD in the UK-ROI cohort. In total, 6582 SNPs within telomere-related genes were present within the discovery GWAS. Comparison of 823 individuals with DKD (cases) with 903 individuals with T1D but no evidence of kidney disease (controls), corrected for age, sex, and duration of diabetes, revealed no genomewide significant SNPs; however, 28 SNPs in 14 genes demonstrated nominal association ( Table 2). The highest-ranked SNP was rs2299694 in ADA (p = 9.77 × 10 −5 , OR = 1.83, 95% CI = 1.35-2.47). All SNPs reaching nominal association underwent in silico replication via meta-analysis with the FinnDiane (n = 2910) and US GoKinD (n = 1595) datasets. Association was supported for two SNPs (Table 2), the most significant being rs2292681 in RNF10 (p = 2.81 × 10 −3 ). Table 2. SNPs in telomere-related genes nominally associated with diabetic kidney disease (DKD) in the discovery cohort, with their association with DKD in the replication cohort also shown. Highlighted in green are those genes significantly associated with DKD in the replication cohort. CHR = chromosome, BP = base pair, A1 = allele 1, OR = odds ratio, SE = standard error, L/U95 = lower/upper 95% confidence intervals.  Analysis for the ESKD phenotype in 247 cases and 1479 controls did not reveal significant SNPs; however, 45 SNPs in 17 unique genes were nominally associated. In silico replication supported nine of these SNPs (Table 3), including SNPs within the top-ranked genes, PIPOX (rs7220474, p = 0.047) and DPP3 (rs2279863, p = 0.028 and rs1671063, p = 0.013). Table 3. SNPs in telomere-related genes nominally associated with end-stage kidney disease (ESKD) in the discovery cohort, with their association with ESKD in the replication cohort also shown. Highlighted in green are those genes significantly associated with ESKD in the replication cohort. CHR = chromosome, BP = base pair, A1 = allele 1, OR = odds ratio, SE = standard error, L/U95 = lower/upper 95% confidence intervals.

Exploring Mendelian Randomisation to Inform Associations between Genetic Telomere Length and DKD
Existing meta-analysis of genetic variants associated with kidney disease in T1D was used as the outcome dataset, considering either DKD or ESKD. As instrumental variables (IVs), 130 variants shown by Codd et al. to be significantly associated with leukocyte telomere length were utilised [48]. Separately, 33 variants used by Park et al. to highlight an association between telomere attrition and CKD were harnessed [47,50]. When utilising either set of IVs, no significant association between telomere length and DKD was identified, with no significant pleiotropy present (Tables 4 and S6). Whilst no significant association was identified between ESKD and telomere length, significant pleiotropy was identified when utilising the Codd et al. IVs (Tables 4 and S6). The MR-PRESSO method confirmed the presence of pleiotropy (global test p = 0.0193) [79]; however, no outlying variants were identified and the causal effect was not significantly distorted.

Significant Differential Methylation of Telomere-Related Genes in DKD and ESKD
Focusing on epigenetics, 496 methylation sites in 212 genes reached epigenome-wide significance (p ≤ 10 −8 ) for DKD (Table S7). cg00445824 in ISYNA1 was the most significantly associated site (p = 9.1 × 10 −24 ), with this site also significant and in the same direction for ESKD. ESKD was associated with 412 methylation sites in 193 unique genes at epigenome-wide significance (Table S8). The most significant was cg19898668 in REM2 (p = 2.2 × 10 −21 ), with this site also significant in the same direction for DKD. Four genes (MAD1L1, TBCD, BANP and PFKB) contained significant differential methylation at more than 10 sites for DKD and ESKD (Table 5). ESKD beta value distributions of top-ranked sites in these genes are shown in Figure 3. Venn diagrams comparing the CpG sites ( Figure S2A) and associated genes ( Figure S2B) show that, whilst 40% of differentially methylated sites overlap between DKD and ESKD, 70% of differentially methylated genes overlap between these phenotypes. Of the top 15 genes presenting differential methylation in both phenotypes (Table 5), MAD1L1, TUBB, HIST1H2AL and TBCA were significantly associated with diabetes-related phenotypes in the UK Biobank (Table 6). MAD1L1 was fine mapped for top-ranking methylation sites ( Figure 4). Diabetes mellitus in pregnancy  Gene ontology analysis of the biological processes enriched in the significantly differentially methylated telomere-related genes, compared to the full telomere-related profile, revealed associations with telomere and chromosomal maintenance. Differential methylation within genes with predicted roles in Wnt signalling was also observed for DKD and ESKD (Tables S9 and S10). Morphogenesis, biosynthetic processing and metabolism regulation were additionally highlighted for DKD alone (Table S9).  Diabetes mellitus in pregnancy  Gene ontology analysis of the biological processes enriched in the significantly ferentially methylated telomere-related genes, compared to the full telomere-related file, revealed associations with telomere and chromosomal maintenance. Differe methylation within genes with predicted roles in Wnt signalling was also observe DKD and ESKD (Tables S9 and S10). Morphogenesis, biosynthetic processing and m olism regulation were additionally highlighted for DKD alone (Table S9). Gene ontology analysis of the biological processes enriched in the significantly differentially methylated telomere-related genes, compared to the full telomere-related profile, revealed associations with telomere and chromosomal maintenance. Differential methylation within genes with predicted roles in Wnt signalling was also observed for DKD and ESKD (Tables S9 and S10). Morphogenesis, biosynthetic processing and metabolism regulation were additionally highlighted for DKD alone (Table S9).

Gene Expression Changes in Telomere-Related Genes Reflected Differential Methylation Patterns
To explore the effects of differential methylation during DKD, methylation data were compared with gene expression data from micro-dissected glomerular and tubulointerstitial tissue [71]. The resultant log fold change (logFC) data from a differential expression analysis in DKD, compared to living donor controls, were correlated with mean delta-beta values for the significantly differentially methylated CpG sites within each gene ( Figure 5). Some genes displayed a change in methylation consistent with their expression pattern during DKD (with fold change greater than 2 in the increasing or decreasing direction); Figure 6 explores these genes further. Interestingly, a number of genes which presented directional change in transcriptomic data, which was concordant with the differential methylation in DKD or ESKD, were associated with hypertension, diabetes, renal, glomerulonephritis or nephritis-related phenotypes in the UK Biobank (Table 7).

Gene Expression Changes in Telomere-Related Genes Reflected Differential Methylation Patterns
To explore the effects of differential methylation during DKD, methylation data were compared with gene expression data from micro-dissected glomerular and tubulointerstitial tissue [71]. The resultant log fold change (logFC) data from a differential expression analysis in DKD, compared to living donor controls, were correlated with mean delta-beta values for the significantly differentially methylated CpG sites within each gene ( Figure  5). Some genes displayed a change in methylation consistent with their expression pattern during DKD (with fold change greater than 2 in the increasing or decreasing direction); Figure 6 explores these genes further. Interestingly, a number of genes which presented directional change in transcriptomic data, which was concordant with the differential methylation in DKD or ESKD, were associated with hypertension, diabetes, renal, glomerulonephritis or nephritis-related phenotypes in the UK Biobank (Table 7).  Only genes shown to be significantly dysregulated (p-adj < 0.05) are displayed. Shown in black are genes of interest where delta-beta methylation values correlate accordingly with expression level changes (for example, increased methylation correlates with decreased expression); shown in red are genes which do not. Labelled are genes with ±1 logFC.

× 10
Essential hypertension values for the CpG sites within these genes of interest, whose differential methylation was significantly associated with DKD or end-stage kidney disease (ESKD) (C). ns denotes genes that were not significantly differentially methylated in DKD/ESKD. Mean ± SD is shown. In graphs A and B, padjusted values from the differential expression analysis are shown.
An additional RNA-sequencing dataset generated from the whole-kidney biopsies of control, early DKD and advanced DKD participants was harnessed [72]. LogFC values for genes of interest in this analysis were generally lower and, therefore, a fold change cut-off of greater than ±1 was utilised ( Figure S3). CAMK1D, ENO2, FANCA, MTX1 and DOK2 displayed consistency between datasets ( Figures S3, S4, 5 and 6). For HSPA6, decreased expression occurred during early DKD, mirroring the first dataset; however, this study revealed a subsequent increase in HSPA6 expression during advanced DKD, which may values for the CpG sites within these genes of interest, whose differential methylation was significantly associated with DKD or end-stage kidney disease (ESKD) (C). ns denotes genes that were not significantly differentially methylated in DKD/ESKD. Mean ± SD is shown. In graphs A and B, p-adjusted values from the differential expression analysis are shown. Table 7. Genes which presented directional change in transcriptomic data which were concordant with the differential methylation in diabetic kidney disease (DKD) or end-stage kidney disease (ESKD) were associated with hypertension, diabetes, renal, glomerulonephritis or nephritis-related phenotypes in the UK Biobank.

Gene
Top An additional RNA-sequencing dataset generated from the whole-kidney biopsies of control, early DKD and advanced DKD participants was harnessed [72]. LogFC values for genes of interest in this analysis were generally lower and, therefore, a fold change cut-off of greater than ±1 was utilised ( Figure S3). CAMK1D, ENO2, FANCA, MTX1 and DOK2 displayed consistency between datasets ( Figures S3 and S4, Figures 5 and 6). For HSPA6, decreased expression occurred during early DKD, mirroring the first dataset; however, this study revealed a subsequent increase in HSPA6 expression during advanced DKD, which may better reflect the decreased methylation associated with ESKD.
When assessing the expression level of telomere-related genes with the highest number of statistically significant methylation sites (Table 5) or those genes with predicted association with Wnt signalling, all genes except WIPI2 showed a significant differential expression in at least one comparison ( Figure S5, Table S11 and Figure S6, Table S12A). Interestingly, C12orf43, TBL1X and TNKS were significantly associated with diabetes or hypertension-related phenotypes in the UK Biobank (Table S12B).

Discussion
In the present study, we explored the association of telomere length with DKD and ESKD in a European T1D population, determining RTL in 1147 participants. Mean RTL was significantly shorter in individuals with DKD, compared to individuals with T1D and no evidence of kidney disease, even after covariate correction. Whilst shorter telomere length has been associated with a significantly increased risk of DKD in diabetes (T1D and T2D combined) [54], mixed associations between telomere length and renal function have been reported. A systematic review by Ameh et al. has investigated associations between telomere length and renal traits in 7,829 individuals from nine studies, two specifically investigating diabetic patients (T1D [53] and T2D [25]) with varying stages of kidney disease [20,25,53]. Telomere attrition was associated with estimated glomerular filtration rate (eGFR) decline and kidney disease progression among people with diabetes [20]. However, longer telomeres were associated with longer kidney disease duration, perhaps influenced by telomere repair in longer surviving CKD patients [20]. Fyhrquist et al. have highlighted that short telomeres were independent predictors of DKD progression in patients with T1D [53]. A recent longitudinal study by Syreeni et al. has revealed that individuals with T1D with the shortest telomeres had lower eGFR, increased albuminuria and more stage 3 CKD [59]. Januszewski et al., however, determined that, whilst telomeres were shorter in patients with T1D versus a control group, RTL did not correlate significantly with renal function [58]. A recent review from our group highlights the recent literature connecting telomere length and DKD [80].
Telomere length does not correlate well with all aspects of renal function, especially after age adjustment [21,22,49,55,81,82]. Sun et al. determined that SNPs in telomere-related genes, rather than telomere length, contributed to primary glomerulonephritis/ESKD susceptibility [49]. Other studies determined that genetic telomere length was not significantly associated with CKD or diabetes, perhaps due to residual biases or limited power [44,47]. These data reflect that telomere length is a weak marker of ageing, displaying substantial inter-individual variation, reflecting differing exposomes [83]. For example, variation in ageing linked to renal function is also strongly influenced by diet and the microbiome, which may confound primary analyses [84]. This highlights the need for multi-omic studies to assess additional factors, such as epigenetics, which can be altered during the life course [5,64,85]. We therefore conducted, to our knowledge, the first investigation of leukocyte telomere length, together with both genetic and epigenetic status of nuclear telomere-related genes, for association with DKD in T1D. We used costeffective methods to determine telomere length and high-density microarrays for efficient and high-throughput examination of telomere-related gene data extracted from~1 million SNPs and~480,000 methylation sites at single-base resolution.
Using existing genotype data [62], SNPs were investigated for an association with qPCR-established telomere length, DKD and ESKD. A total of 17% of the top 36 SNPs associated with telomere length were present in WRN. Variants in this gene are the primary cause of Werner syndrome, a disorder of accelerated ageing, with WRN knockout cell lines presenting with accelerated telomere attrition [86][87][88]. In vitro studies have demonstrated that the WRN helicase interacts with TERF2, a member of the shelterin complex essential for telomere maintenance [87,89]. TERF1 is an additional protein implicated in WRN function, with a variant near TERF1-interacting nuclear factor 2 (TINF2) associated with telomere length and T2D, and is shown to drive the association between telomere length and CKD in T2D [51,52,90]. Together, these studies highlight the potential for WRN functioning and TERF1-or TERF2-interacting proteins to influence the risk of DKD. It is important to note, however, that no SNP reached genome-wide significance when assessing RTL, and genes such as LMNA, implicated in telomere stability and CKD [91,92], were not identified, highlighting the need for even larger genetic datasets of this type to increase the power and uncover novel gene-function interactions.
The nominal association with DKD was supported during replication analysis for two SNPs, the most significant being rs2292681 in RNF10. Liu et al. recently showed via a multi-ancestry meta-analysis of 1.5 million individuals that rs3817301 in RNF10 was significantly associated with eGFR [93]. For ESKD, variants within seven genes were supported during replication analysis, with variants within SYK and PIPOX previously associated with T2D or HbA1c [94,95]. PIPOX expression was shown to significantly decrease during DKD [96], encoding a protein involved in maintaining oxidative stress balance [96], a process disrupted during DKD [97] and a key factor in modulating telomere shortening [98].
Utilising Mendelian randomisation, we assessed the causal association between genetically determined telomere shortening and DKD or ESKD, harnessing two sets of IVs previously used to explore kidney disease in European populations [47,48,50]. Park et al. determined that telomere shortening was significantly associated with a higher CKD risk (OR 1.20, p < 0.001), with successful replication in the UK Biobank [47,50]. Li et al. utilised the same genetic instrument in a UK Biobank cohort to reveal that geneticallydetermined telomere attrition did not affect the risk of diseases such as diabetes or CKD [47]; however, Park et al. highlighted that this null result may be due to the low number of self-reported or ICD-confirmed diagnosis for CKD and instead opted for serum cystatin C and creatinine level classification [50]. Codd et al. identified IVs of leukocyte telomere length in a UK Biobank cohort and, whilst these authors determined that direct measurements of telomere length were significantly associated with CKD within the UK Biobank (Hazard ratio = 0.889, p value = 9.45 × 10 −17 ), genetically determined telomere length was not (OR = 1.02, p = 0.71), with CKD classified based on self-reported, ICD and procedure codes [48]. Harnessing these genetic instruments for analysis within our cohort revealed no significant associations between telomere shortening and DKD or ESKD. Additional studies have explored the effect of genetically determined telomere length and kidney disease. Gurung et al. demonstrated that genetically determined telomere attrition was associated with increased DKD in East Asian T2D patients [52]. Taub et al. highlighted that rs1008438 in HSPA1A gene was significantly associated with DKD risk in T1D [46]. These results highlight the complexity of the association between telomere length and kidney disease, with results in both genetic-based and observational studies dependent on the measure of kidney function used [22]. It is important, therefore, to take a wider approach, assessing both the observational and genomic nature of telomere length to gain a fuller understanding of kidney disease pathogenesis.
Via a multi-omics approach, we broadened our study of telomere-related genes in DKD. Multiple statistically significant methylation sites were identified in MAD1L1, encoding a protein involved in the mitotic spindle assembly, chromosome alignment, cell cycle control and tumour suppression. MAD1L1 inhibits TERT transcription via epigenetic modification of histones [99]. TERT encodes the catalytic subunit of telomerase and plays a crucial role in telomere maintenance and senescence [100]. Methylation of TERT and MAD1L1 in DKD may alter telomerase activity, potentially reducing the regenerative ability of renal cells [101,102]. Via TERT, telomerase can also modulate Wnt/β-catenin signalling [103]. This pathway is important for podocyte proliferation, the epithelial cells involved in maintaining normal kidney filtration (Figure 7). Interestingly, gene ontology analysis revealed that Wnt signalling was an enriched biological process within the telomererelated genes which were significantly differentially methylated in DKD or ESKD. Altering DNA methylation, even via dietary interventions [84,104], may prove a novel therapeutic opportunity for DKD [105], with this analysis identifying targets for future study (TERT, VAX2, C12orf43, TBL1X, TNKS, BCL7B, ZEB2 and AKT1). . Via TERT, telomerase can also modulate Wnt/β-catenin signalling, with this mechanism controlling proliferation independent of TERT catalytic activity or the telomerase RNA component (TERC) [103,106,107]. Wnt signalling influences on kidney development and functioning, with potential for this pathway to be targeted for therapeutic interventions for DKD [108][109][110][111][112][113][114][115][116][117].
A limitation of this study was the analysis of telomere length and epigenetics at a single time point. A longitudinal study, as carried out by Syreeni et al. in a smaller T1D cohort [59], would provide insights into disease progression and (epi)genomic changes over time. In addition, incorporating assessment of individual exposome factors, such as lifestyle, diet, environment and socioeconomic position, as well transgenerational expo- Figure 7. Via TERT, telomerase can also modulate Wnt/β-catenin signalling, with this mechanism controlling proliferation independent of TERT catalytic activity or the telomerase RNA component (TERC) [103,106,107]. Wnt signalling influences on kidney development and functioning, with potential for this pathway to be targeted for therapeutic interventions for DKD [108][109][110][111][112][113][114][115][116][117].
A limitation of this study was the analysis of telomere length and epigenetics at a single time point. A longitudinal study, as carried out by Syreeni et al. in a smaller T1D cohort [59], would provide insights into disease progression and (epi)genomic changes over time. In addition, incorporating assessment of individual exposome factors, such as lifestyle, diet, environment and socioeconomic position, as well transgenerational exposome stressors, may also be of merit, due to their impact on accelerated ageing, kidney disease and the epigenome [83].
Employing previously published RNA-sequencing datasets generated using human kidney tissue biopsies [71,72], differential methylation of telomere-related genes during DKD and ESKD was correlated with altered gene expression. Differentially methylated genes with predicted roles in Wnt signalling were significantly differentially expressed during DKD; however, the extent of this differential expression was limited. Differentially methylated genes showing a concordant change in gene expression during DKD were identified. Genetic variants within BICD1 were nominally associated with DKD and ESKD in the discovery and replication GWAS, with rs7900065 within BAG3 nominally associated with ESKD in both datasets. A BICD1 genetic variant was previously associated with longitudinal DNA methylation changes in an older population [78]. The present study also highlighted a telomere-related gene (SH3GL2), which contained genetic variants nominally associated telomere length and genetic variants within 150 kilobases previously associated with DNA methylation or accelerated epigenetic ageing. Interestingly, SH3GL2 was significantly differentially methylated in ESKD in this study and significantly dysregulated in both the Fan et al. and Levin et al. RNA-sequencing datasets, with decreased expression observed alongside (on average) decreased methylation. These results highlight the potential interaction between genetic variation in telomere-related genes, methylation changes and gene expression, which may impact kidney health during T1D, identifying potential diagnostic and therapeutic targets.
Whilst a number of genes with a high number of differentially methylated sites, such as MAD1L1, presented significant differential expression, the extent of this differential expression was limited. Therefore, specifically methylated sites, rather than an accumulation of differentially methylated sites, may influence gene expression during DKD. Future work is needed to explore the effects of differential methylation in DKD and ESKD, especially for genes where extensive transcriptional change was not observed or where gene expression patterns were not concordant with methylation changes.

Conclusions
In summary, taking a multi-omic approach, telomere-related genes nominally associated with telomere length, DKD, or ESKD were identified, with epigenetic analysis revealing approximately 230 genes differentially methylated in DKD and/or ESKD. These differentially methylated genes were enriched for Wnt signalling functions. Harnessing previously published transcriptomic datasets, differential methylation was correlated with gene expression in DKD, highlighting potential targets where epigenetic dysregulation may result in altered gene expression and influence disease. This study identified potential targets where epigenetic regulation of telomere function may have functional consequence on DKD, useful as potential diagnostic and therapeutic targets for intervention.  . Figure S1. UK ROI GWAS to determine statistically significant genetic associations for relative telomere length. (A) QQ plot highlighted limited deviation of the observed p values from the null hypothesis. (B) Manhattan plot with those SNPs reaching the −log10(p) > 5 cut-off annotated. (C) Zoomed Manhattan plot displaying those SNPs reaching the p < 0.005 cut-off. Shown in red are SNPs within the 376 telomere-related genes. The leading SNP within each telomere-related gene is annotated. Figure S2. Venn diagrams to highlight the number of differentially methylated CpG sites overlapping between diabetic kidney disease (DKD) and end-stage kidney disease (ESKD) (A) or the overlapping genes associated with these differentially methylated CpG sites (B). Figure S3. Mean delta-beta (∆β) values per gene were determined for the differentially methylated CpG sites significantly associated with diabetic kidney disease (DKD) (A,C,E) or end-stage kidney disease (ESKD) (B,D,F). Mean ∆β values were regressed on log fold change (logFC) values obtained during the differential expression analysis of early DKD versus control (A,B), advanced DKD versus control (C,D) or advanced DKD versus early DKD (E,F) RNA-sequencing data [72]. Only genes shown to be significantly dysregulated (p-adj < 0.05) are displayed. Shown in black are genes of interest where ∆β values correlate accordingly with expression level changes (for example, increased methylation correlates with decreased expression); shown in red are genes which do not. Labelled are genes with ±0.1 logFC. Figure S4. Expression level (normalised log2 read count) for control, early diabetic kidney disease (DKD) or advanced DKD (Fan et al.) [72] for genes of interest (A). P-adjusted values from the differential expression analysis are shown in star format: * p < 0.05, ** p < 0.01, *** p < 0.001 (nonsignificant differences are not displayed). Below are the delta-beta (∆β) values for the same genes which were significantly differentially methylated in DKD or end-stage kidney disease (ESKD) phenotypes (B). Note that, for ARK1B10 and DOK2, delta-beta values are not shown for the DKD phenotype as changes in methylation within these genes was not significantly associated with DKD. For AHCY, ARRK1 and TWF2, delta-beta values are not shown for the ESKD phenotype are not shown as changes in methylation within these genes were not significantly associated with ESKD. ns denotes genes that were not significantly differentially methylated in DKD/ESKD. Mean ± SD is shown. Figure S5. Expression level (normalised log2 read count) for control and diabetic kidney disease (DKD) tubulointerstitial (A) or glomerular (B) tissue (Levin et al.) [71], for genes with the highest number of significantly differentially methylated sites in DKD or end-stage kidney disease (ESKD) phenotypes. Also shown are expression levels (normalised log2 read count) for control, early DKD or advanced DKD in an alternative dataset (Fan et al.) [72] (C). p-adjusted values from the differential expression analysis are shown. Delta-beta (∆β) values for the same genes which were significantly differentially methylated in DKD or ESKD phenotypes (D). Mean ± SD is shown. Figure S6. Expression level (normalised log2 read count) for control and diabetic kidney disease (DKD) tubulointerstitial (A) or glomerular (B) tissue (Levin et al.) for genes having a predicted association with Wnt signalling via gene ontology analysis. Also shown are expression levels (normalised log2 read count) for control, early DKD or advanced DKD in an alternative dataset (Fan et al.). (C) p-adjusted values from the differential expression analysis are shown. Delta-beta values for the same genes which were significantly differentially methylated in DKD or end-stage kidney disease (ESKD) phenotypes (D). Mean ± SD is shown. Note that, for TERT and VAX2, expression values were not determined in either of the RNA-sequencing datasets analysed. Table S1. Genes with a known role in telomere function. SNPs for these relevant genes were downloaded from the Ensembl genome browser using Ensembl genes 74 database on the Homo sapiens dataset. Shown here are the gene start and end positions (base pairs (bp)) in the GRCh37.p13 assembly. Table S2. Primer sequences for fine mapping of selected topranked CpG. Table S3. Genetic variants in the UK-ROI population nominally associated with telomere length. NA = no gene assigned, CHR = chromosome, BP = base pair, A1 = allele 1, OD = odds ratio, SE = standard error, L95 = beta lower 95% confidence limit, U95 = beta upper 95% confidence limit, STAT = test statistic, p = p-value. Table S4. SNPs in telomere-related genes (Table S1) demonstrating nominal association (p < 0.005) with telomere length. Also shown are the number (or percentage) of the top 20 phenotypes for each variant, identified via PheWEB, associated with hypertension, diabetes, renal, glomerulonephritis or nephritis-related phenotypes. CHR = chromosome, BP = base pair, A1 = allele 1, SE = standard error, L95 = lower 95% confidence limit, U95 = upper 95% confidence limit, STAT = test statistic, p = p-value. Table S5. SNPs in telomere-related genes (Table S1) demonstrating nominal association (p < 0.005) with telomere length (Table S4) and the SNP regions and hits for renal, diabetes and cardiovascular phenotypes on the Common Metabolic Disease (CMD) Portal. Also shown are hits within these regions in the GWAS Catalogue for telomere biology, longevity or methylation phenotypes. NA = not applicable, CHR = chromosome, BP = base pair. Table S6. Heterogeneity statistics for MR analysis with the replication meta-analysis data demonstrating the effect of telomere attrition on the risk of diabetic kidney disease or end-stage kidney disease in T1D, using the Codd et al. instrumental variables (IVs) or the Park et al. IVs [47,48,50]. IVW = inverse variance weighted. DKD = diabetic kidney disease. ESKD = end-stage kidney disease. Table S7. The 496 CpG sites in 212 unique genes with a known role in telomere function demonstrating genomewide significance with diabetic kidney disease in the UK-ROI population. Arranged in order of significance. Table S8. The 412 CpG sites in 193 unique genes with a known role in telomere function demonstrating genome-wide significance with end-stage kidney disease in the UK-ROI population.
Arranged in order of significance. Table S9. Clustered significantly enriched biological process gene ontology terms for the significantly differentially methylated telomere-related genes associated with diabetic kidney disease. Gene frequency is shown as a percentage and fraction (number of significant genes/number of annotated genes). Table S10. Clustered significantly enriched biological process gene ontology terms for the significantly differentially methylated telomere-related genes associated with end-stage kidney disease. Gene frequency is shown as a percentage and fraction (number of significant genes/number of annotated genes). Table S11. Differential expression results for the genes presenting the highest number of differentially methylated CpG sites associated with diabetic kidney disease (DKD) or end-stage kidney disease. Log fold change (logFC) and adjusted p-values (adj p-value) are shown. Significant differential expression is highlighted using an asterisk (*). Table S12. (A) Differential expression results for the genes containing differentially methylated CpG sites associated with diabetic kidney disease (DKD) or end-stage kidney disease which are predicted to associate with Wnt signalling. Log fold change (logFC) and adjusted p-values (adj p-value) are shown. Significant differential expression is highlighted using an asterisk (*). (B) Genes from Table S12(A) which present associations with hypertension, diabetes, renal, glomerulonephritis or nephritis-related phenotypes, as determined via PheWEB. Shown is the top p-value per association per gene. Institutional Review Board Statement: Ethical review was obtained from district ethics committees for individual studies, as described previously for these individual cohorts [62,119,121,122].
Informed Consent Statement: Informed consent was obtained from all individuals involved in the study.

Conflicts of Interest:
P.-H.G. reports receiving lecture honorariums from Astellas, AstraZeneca, Bayer, Boehringer Ingelheim, Eli Lilly, Elo Water, Medscape, MSD, Mundipharma, Novo Nordisk, PeerVoice, Sanofi, Sciarc, and being an advisory board member of AbbVie, Astellas, AstraZeneca, Bayer, Boehringer Ingelheim, Eli Lilly, Medscape, MSD, Mundipharma, Nestlè, Novo Nordisk, and Sanofi. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript; or in the decision to publish the results.