Early DNA methylation changes in children developing beta cell autoimmunity at a young age

Aims/hypothesis Type 1 diabetes is a chronic autoimmune disease of complex aetiology, including a potential role for epigenetic regulation. Previous epigenomic studies focused mainly on clinically diagnosed individuals. The aim of the study was to assess early DNA methylation changes associated with type 1 diabetes already before the diagnosis or even before the appearance of autoantibodies. Methods Reduced representation bisulphite sequencing (RRBS) was applied to study DNA methylation in purified CD4+ T cell, CD8+ T cell and CD4−CD8− cell fractions of 226 peripheral blood mononuclear cell samples longitudinally collected from seven type 1 diabetes-specific autoantibody-positive individuals and control individuals matched for age, sex, HLA risk and place of birth. We also explored correlations between DNA methylation and gene expression using RNA sequencing data from the same samples. Technical validation of RRBS results was performed using pyrosequencing. Results We identified 79, 56 and 45 differentially methylated regions in CD4+ T cells, CD8+ T cells and CD4−CD8− cell fractions, respectively, between type 1 diabetes-specific autoantibody-positive individuals and control participants. The analysis of pre-seroconversion samples identified DNA methylation signatures at the very early stage of disease, including differential methylation at the promoter of IRF5 in CD4+ T cells. Further, we validated RRBS results using pyrosequencing at the following CpG sites: chr19:18118304 in the promoter of ARRDC2; chr21:47307815 in the intron of PCBP3; and chr14:81128398 in the intergenic region near TRAF3 in CD4+ T cells. Conclusions/interpretation These preliminary results provide novel insights into cell type-specific differential epigenetic regulation of genes, which may contribute to type 1 diabetes pathogenesis at the very early stage of disease development. Should these findings be validated, they may serve as a potential signature useful for disease prediction and management. Graphical abstract Supplementary Information The online version contains peer-reviewed but unedited supplementary material available at 10.1007/s00125-022-05657-x.


Introduction
Type 1 diabetes is a complex autoimmune disease with a strong genetic component. Genome-wide association studies (GWAS) have identified more than 60 genetic loci associated with the risk of type 1 diabetes [1][2][3]. However, genetic variation alone cannot explain the conspicuous rise in the disease incidence during the past decades [4]. Several environmental factors, including viral infections, diet, toxins and other determinants, have been implicated [5].
Environmental exposures may result in epigenetic modifications of DNA, chromatin and histone proteins that change the level of gene expression. The importance of epigenetic gene regulation in complex disease phenotypes has been highlighted for several autoimmune disorders, such as rheumatoid arthritis (RA) and inflammatory bowel disease (IBD) [6,7].

Coverage filtering
Removed CpGs with coverage >99.9th percentile Merged technical replicates Required minimum coverage 10 in at least 2 matching time points in at least 5 (out of 7) case-control pairs
A few studies identified a link between DNA methylation in circulating immune cells and type 1 diabetes [9][10][11][12][13]. However, most of those studies focused on individuals who had already developed clinical type 1 diabetes. A recent study examined DNA methylation changes in peripheral blood samples longitudinally collected before the diagnosis of type 1 diabetes from children enrolled in the prospective Diabetes Autoimmunity Study in the Young (DAISY) cohort [14]. However, the study was limited by the use of a 450 K array for methylation profiling that restricts the analysis to approximately 450,000 targeted CpG sites.
The detection of type 1 diabetes-associated autoantibodies targeting insulin (IAA), glutamic acid decarboxylase (GADA), tyrosine phosphatase-like protein (islet antigen-2, IA-2A) and zinc transporter-8 (ZnT8A) is currently the primary method to predict the development of type 1 diabetes. Children who develop multiple autoantibodies are very likely to progress to type 1 diabetes [15]. However, the appearance of autoantibodies is an indication of an already ongoing autoimmune reaction. Therefore, discovery of methylation changes at the early stages of disease onset, before seroconversion, could provide new insights into molecular mechanisms leading to type 1 diabetes.
The aim of this study was to assess early DNA methylation changes associated with type 1 diabetes in longitudinally collected samples from children who later developed islet autoimmunity at a young age. We used the reduced representation bisulphite sequencing (RRBS) approach instead of array-based methods to achieve broader coverage including sites that are not included in the arrays. Further, we purified the CD4 + T cell, CD8 + T cell and CD4 − CD8 − cell fractions of peripheral blood mononuclear cells (PBMCs) to determine cell type-specific methylation differences that have not been reported earlier.
Besides identifying type 1 diabetes-associated DNA methylation changes, we investigated correlations between DNA methylation and gene expression, exploiting our earlier RNA sequencing (RNA-seq) study of the same samples [16]. The identification of early DNA methylation changes associated with type 1 diabetes will provide novel insights into the pathogenesis of the disease, and may serve as a potential signature for disease prediction and management.

Methods
Study cohort The cohort analysed here using RRBS is described in an earlier report on transcriptomics analysis of the same individuals [16]. DNA samples were longitudinally collected from seven case-control pairs of the Pathogenesis of Type 1 Diabetes -Testing the Hygiene Hypothesis (DIABIMMUNE) cohort (Fig. 1a, electronic supplementary material [ESM] Table 1). In total, we analysed 226 samples: 73, 77 and 76 for CD4 + T cells, CD8 + T cells and CD4 − CD8 − cell fractions, respectively. The study protocols were approved by the ethics committees of the participating hospitals, and the parents gave written informed consent. HLA genotyping was described earlier [16]. Autoantibodies IAA, GADA, IA-2A and ZnT8A were measured from serum with specific radiobinding assays. The cutoff values were based on the 99th percentiles in children without diabetes, which were 2.80 relative units (RU) for IAA, 5.36 RU for GADA, 0.78 RU for IA-2A and 0.61 RU for ZnT8A. A sample was considered seropositive when any of the autoantibodies exceeded the threshold. The case individuals became positive for at least two type 1 diabetes-specific autoantibodies at the age of 1-2 years, whereas the control individuals remained autoantibody-negative throughout the follow-up period. The case-control pairs were matched by sex, place of birth, age and HLA risk class. The age was matched ±2 months except for the 36 months samples, where the samples were collected within 4 months. The HLA risk classes were defined as described earlier [17]. A detailed description of the study participants and their antibody measurements is presented in ESM Table 2.

RRBS library preparation
The library preparation from 100 ng of genomic DNA was carried out according to a protocol adapted from a gel-free multiplexed RRBS method [18]. Subsequent steps were performed as described earlier [19].  DMCs. Methylation profiles of top-ranked DMCs between cases and controls identified in CD4 + T cells (a-d), CD8 + T cells (e, f) and CD4 − CD8 − fractions (g, h). The plots are coloured according to the pairwise methylation difference, as specified in the colour key. Blue and red indicate hypo-and hypermethylation of the CpG site in cases, compared with controls, respectively. Black represents missing values. The pairs of individuals are not in any particular order. A black dot in the boxes indicates a sample collected after seroconversion. m, months 2 × 100-bp read length was used with a 6-bp index run. Technical quality of the HiSeq 2500 run was good, and the cluster amount was as expected. More than 76% of all bases above Illumina quality score Q30 was required. The sequencing runs of 37 samples were repeated due to low raw read counts (less than 30 million reads per sample). The median yield of the other samples was 47 million reads (read 1 + read 2).
RRBS data analysis RRBS data analysis was performed as we described previously [20]. The data analysis workflow is summarised in Fig. 1b. More detailed description is also provided in the ESM Methods. Most of the data analysis was done with R versions 3.6.1 and 4.0.4 [21]. The heatmaps in Fig. 2 were produced with the R package hamlet [22]. Annotation of differentially methylated CpG sites (DMCs) to genomic regions was carried out through R package genomation version 1.16.0 [23] using Genome Reference Consortium Human Build 37 (GRCh37/hg19).
Expression quantitative trait methylation analysis A list of human whole blood cis-expression quantitative trait methylation (cis-eQTM, false discovery rate-corrected p value [FDR] <0.1) effects was downloaded from the BIOS QTL database [24] (accessed 1 April 2020), which is based on 3841 Dutch peripheral blood samples.
Methylation-expression correlation analysis We utilised our published RNA-seq data from the same samples to calculate methylation-expression correlations. Spearman correlations across all samples within each cell fraction were calculated between each DMC and all genes with transcription start sites (TSSs) within a window of 250 kb in both directions from the CpG site's genomic location. The Ensembl identifiers and gene names are from Ensembl versions 70 [25] and 75 [26].
Pathway enrichment analysis Gene Ontology (GO) terms were used in the pathway enrichment analysis. Only DMCs with adjusted p value (Benjamini-Hochberg-corrected p value) <0.1 and their corresponding nearest and correlating genes were used for the analysis. Pathways were considered to be significantly enriched at FDR <0.05 (Fisher's exact test).
Pyrosequencing PyroMark assay design 2.0 software (Qiagen) was used to design the assays. Target-specific primers designed with the software include a primer pair with site-specific biotinylation for target amplification and a pyrosequencing primer. Sample preparations started with 200 ng of DNA from three selected case-control pairs. Samples were first sodium bisulphite-treated with EZ DNA Methylation-Gold Kit (Zymo Research, Irvine, CA, USA, cat. No. D5006). Selected targets were then amplified by PyroMark PCR Kit (Qiagen, cat. No.

Results
RRBS analysis of different cellular fractions identified methylation differences in children at risk for type 1 diabetes When including all the longitudinal samples into the differential methylation analysis model, we discovered 225, 114 and 87 DMCs in the CD4 + T cell, CD8 + T cell and CD4 − CD8 − cell fractions, respectively, at FDR <0.1 and absolute coverage-corrected mean methylation difference >0.1 (Fig. 2, Table 1, ESM Table 3-5).
In the CD4 + T cells, the DMC located in the intergenic region between TRAF3 and RCOR1 was hypomethylated in cases. Furthermore, three highly significant DMCs were found at the promoters of protein-coding genes: RGPD8, FOXR1 and ARRDC2. A variant in ARRDC2 has been genetically associated with early-onset Crohn's disease [27]. In CD8 + T cells, the most significant DMC was located in the intergenic region between the pseudogene LOC390705 and a non-coding RNA, LOC113002582. In the CD4 − CD8 − cell fraction, highly significant differences in methylation were observed near proteincoding genes WDR7 and EXOC5 (Table 1).
Further, to gain insights into biological functions of the genes near DMCs, we performed GO pathway enrichment analysis (Fisher's exact test) on the sets of the nearest genes for each cell type. For CD8 + T cells, two GO terms were enriched: peptidyl-tyrosine phosphorylation (FDR = 0.035) and cell differentiation (FDR = 0.045). No GO terms were significantly enriched among the nearest genes for CD4 + T cells or CD4 − CD8 − cell fractions.
Most of the DMRs were cell type-specific Although a single CpG site might affect gene expression, often adjacent CpG sites may act together in gene expression regulation [28]. Such DMCs can be combined into a differentially methylated region (DMR). Therefore, we combined individual DMCs within 2 kb with consistent methylation patterns into DMRs. We identified 79, 56 and 45 DMRs between the cases and controls in the CD4 + T cell, CD8 + T cell and CD4 − CD8 − cell fractions, respectively (Fig. 3a, ESM Table 6 -8). While the majority of DMRs were cell type-specific, three DMRs near the FOXR1, RGPD8 and LOC100128946 (also known as LINC01310) genes were hypermethylated and one DMR near the LOC390705 pseudogene was hypomethylated in cases in all three cell fractions. These results suggest that differences in these three DMRs are conserved across cell types. It is also possible that DNA methylation at these DMRs regulates the expression of genes with a similar functional role in all three subsets. fractions constituted up to 30%. In total, 9-11% of DMRs were located in exons in all three cell fractions. The GeneHancer database was utilised to explore DMRs in the enhancer regions. Ten, four and seven DMRs were located in enhancer regions in CD4 + T cells, CD8 + T cells and the CD4 − CD8 − cell fractions, respectively (ESM Tables 6-8). The DMR near CBFA2T3 is a part of an enhancer that interacts with its promoter based on the GeneHancer database and may regulate its expression (Fig. 4a). CBFA2T3 encodes a transcriptional repressor and plays a role in T cell development. A variant of this gene has been associated with autoimmune vitiligo [29]. We also found DMRs near extended promoter regions of several genes including TOX and IL32 (Fig. 4b, c, ESM Tables 6-8) in the CD4 + and CD8 + T cell fractions, respectively. Interestingly, of the 47, 33 and 17 stand-alone DMCs, which were not part of DMRs, 12, ten and four were in the promoter or enhancer regions in CD4 + T cells, CD8 + T cells and the CD4 − CD8 − cell fractions, respectively (ESM Table 3 -5).
High correlation between methylation and gene expression highlighted candidate genes of interest The availability of RNA-seq data from the purified cell fractions of the very same samples allowed us to assess a possible functional impact of the discovered DMCs on gene expression. We examined correlations between DNA methylation and gene expression, utilising our earlier published RNA-seq data [16] (ESM Tables 3-5, 9 -11). We focused on DMCs that had a high Spearman rank correlation coefficient (>|0.5|) with gene expression (Table 2, Fig. 5). For example, in CD4 + T cells, we found two DMCs on the intron of DGKQ positively correlated with the expression of the gene (Table 2, Fig. 5a). Next, four DMCs on the intron of PCBP3, which is a probable enhancer region, showed an inverse correlation with its expression (Fig. 5b, Table 2). A DMC at the exon of LOC642852 (also known as LINC00205) showed an inverse correlation with the expression of the POFUT2 gene, which is approximately 30 kb away from the DMC. POFUT2 encodes a fucosyltransferase (FUT), responsible for protein O-fucosylation, a post-translational modification process. Interestingly, a genetic risk variant for type 1 diabetes at FUT2, another FUT of a closely related family, was reported previously [30][31][32].
In CD8 + T cells, there were only two DMCs with a Spearman rank correlation >|0.5|. The DMC on chromosome 16 in the intron of LOC100134368 had an inverse correlation with expression of TMEM8A (also known as PGAP6), which is approximately 2 kb away from the DMC (Fig. 5c, Table 2). Another DMC in the exon of ADAM8 correlated inversely with the expression of this gene (Fig. 5d, Table 2).
Further, we analysed CpG sites in our study for already published cis-eQTM associations, using data from a Dutch biobank of cis-eQTM effects based on 3841 whole blood samples [24] (ESM Tables 3-5, [9][10][11]. The expression quantitative trait methylation (eQTM) effects were discovered with the 450 K DNA methylation array, which covers approximately 5% of CpG sites identified in our study. Among those available, we   Tables 3-5, [9][10][11]. The observed direction of expression-methylation correlation in our data was in line with the eQTM effect at CpG sites near PRDM8 and FBXO41.
Pre-seroconversion analysis revealed DNA methylation differences at a very early stage of type 1 diabetes In addition to the differential methylation analysis of all longitudinal samples, we analysed DNA methylation in the samples taken before seroconversion to identify changes that occur prior to detection of islet autoantibodies. We refer to this analysis as pre-seroconversion analysis. Unless specified, the differential methylation analysis refers to the analysis where we take all the longitudinal time points into consideration. We identified 249, 144 and 143 DMCs between cases and controls in the CD4 + T cell, CD8 + T cell and CD4 − CD8 − cell fractions, respectively (FDR <0.1 and coverage-corrected mean methylation difference >0.1) (ESM Tables 9-11). The numbers of DMRs between cases and controls discovered before seroconversion were 97, 68 and 63 in the CD4 + T cell, CD8 + T cell and CD4 − CD8 − cell fractions, respectively (ESM Tables 12-14). Importantly, the overlay of DMRs identified in all longitudinal samples and DMRs discovered prior to the appearance of autoantibodies (Fig. 6a) revealed that the majority of DMRs were pre-seroconversion-specific. We found 58, 52 and 48 pre-seroconversion-specific DMRs in CD4 + T cell, In CD4 + T cells, the most significant DMC was found at the promoter of ARRDC2 gene and was hypomethylated in cases vs controls (Table 3). This DMC was also one of the   (Table 1). In CD8 + T cells, the most significant DMC was observed at the intron of the PCBP3 gene. This CpG site is located in the plausible enhancer region (GeneHancer database) and may regulate the expression of PCBP3. Furthermore, methylation-expression analysis revealed a high correlation between methylation at PCBP3 and expression of the gene (Table 4).
Notably, in CD4 + T cells we identified hypermethylation in case compared with control participants at the promoter of the transcription factor IRF5 and DNAJB12 (Fig. 6b, c, ESM Table 12).
Several DMCs were validated using pyrosequencing analysis in a subset of samples The RRBS results were validated by bisulphite pyrosequencing of the CD4 + and CD8 + T cells of three case-control pairs from the same cohort. We selected nine DMCs in CD4 + T cells and one DMC in CD8 + T cells based on statistical significance and the consistency of pairwise methylation difference profiles.
The direction of methylation difference was concordant between RRBS and bisulphite pyrosequencing results for eight out of ten selected targets (Table 5). High positive correlations between RRBS and pyrosequencing results were observed for CpG sites near ARRDC2 (chr19:18118304), PCBP3 (chr21:47307815) and TRAF3 (chr14:81128398) and also at the intergenic CpG site near PRDM8 (chr4:81128398) with a known eQTM effect on the gene expression. In the CD8 + T cells, we validated the DMC near IL32 (chr16:3116115), which was hypomethylated on the promoter region in type 1 diabetes, compared with control participants.

Discussion
In this study, we identified DNA methylation changes in children developing type 1 diabetes, compared with their matched control participants, before diabetes diagnosis and even before the appearance of autoantibodies associated with type 1 diabetes. Analysing fractionated PBMC samples allowed us to discover immune cell subset-specific DNA methylation changes. Importantly, DNA methylation changes were observed in genes associated with type 1 diabetes, including IL32, TRAF3 and DGKQ, and in novel candidate genes, which have not been linked to the disease, such as ARRDC2 and PCBP3. Using pyrosequencing we validated the direction of methylation difference identified by RRBS for eight out of ten selected targets. Since the number of case individuals in this prospective study cohort was limited (n = 7), we selected matched control individuals to control for effects of confounding factors, such as sex and place of birth. The risk of error was further minimised by analysing several time points for each individual. Even with these measures taken, the number of individuals is an important limitation of this study, and the results need to be validated in a larger cohort.
A strength of our study is the correlation of DNA methylation with gene expression from the purified cell type of the very same samples from our earlier report [16], which enabled us to highlight candidate genes whose expression might be influenced by differential methylation. For instance, an intron of PCBP3 had a locus hypomethylated in cases in CD4 + T cells across all time points and in CD8 + T cells before seroconversion, and it had an inverse correlation with the gene expression. RNA binding proteins of the poly(rC) binding protein (PCBP) family are important in regulating gene expression at the post-transcriptional level. Further, another locus at the promoter of IL32 was also hypomethylated in cases compared with controls in the analysis of all longitudinal samples in CD8 + T cells. Methylation-expression analysis showed a weak inverse correlation between IL32 mRNA expression and promoter methylation, suggesting that epigenetic changes at the promoter of IL32 lead to a higher expression of this cytokine. Interestingly, IL-32 was significantly upregulated at the mRNA level in cases compared with controls in this cohort. IL-32, a proinflammatory cytokine, The first six columns contain information about selected targets. In column 7 we indicate whether the methylation difference between the cases and controls was concordant (to the same direction) in the pyrosequencing results compared with RRBS results. Column 8 contains the Pearson correlation coefficient between pyrosequencing and RRBS results (methylation proportions of the six samples that were pyrosequenced). Columns 9 and 10 contain raw and adjusted p values (FDR), respectively was shown to be upregulated at either mRNA or protein level in several autoimmune disorders [33]; however, the precise role of this cytokine remains to be elucidated. Next, we observed hypermethylation in cases vs controls at a locus near DGKQ, which correlated positively with the expression of this gene in CD4 + T cells and CD4 − CD8 − cell fractions. This finding suggests that the differential methylation at DGKQ may result in its higher expression in type 1 diabetes cases. This observation is in line with the results of the transcriptomic analysis of CD4 + T cells from the same cohort [16]. Diacylglycerol kinase theta (DGKQ) is a member of a family of diacylglycerol (DAG) kinases, which have critical roles in T cell receptor (TCR) signalling [34]. The kinase inversely regulate T cell activation, by terminating DAG-mediating signalling [35]. Interestingly, a recent integrative analysis of genetic association, functional genomics data and protein-protein networks identified DGKQ as a potential novel type 1 diabetes therapeutic gene target along with other novel and known targets, e.g., IL2RA, IL6ST, IL6R and TYK2 [3].
The type 1 interferon transcriptional signature is associated with the development of type 1 diabetes [36,37]. In CD4 + T cells of pre-seroconversion samples, we observed hypermethylation at the promoter of IRF5 (Fig. 6b, ESM Table 12). This DMR is located in the vicinity (less than 7 kb) of multiple SNPs that were previously associated with RA [38,39], systemic lupus erythematosus (SLE) [40,41] and IBD [42]. This finding suggests that the transcription factor may have a role at the very early stages of islet autoimmunity. Further, in CD8 + T cells, we found hypermethylation at an exon of IRF1 gene (ESM Table 7). Interferon regulatory factor 1 (IRF1) is important in the host immune response to pathogens and CD8 + T cell maturation [43]. Further, we found hypomethylated DMRs at TRAF3 in CD4 + T cells and CD4 − CD8 − cell fractions. TRAF3 is also critical in regulation of type 1 interferon production and plays a role in antiviral immune response [44]. In CD4 + T cells, hypomethylation at TRAF3 had moderate positive correlation with the gene expression (Spearman ρ = 0.334), which might result in lower gene expression in cases vs controls. TNF receptor-associated factor 3 (TRAF3) was shown to be recruited in the TCR/CD28 signalling complex in murine CD4 + T cells and critical for T cell-mediated immunity to infection [45].
Thymocyte selection-associated high mobility group box protein (TOX) is a key transcriptional regulator indispensable for CD4 + T cell development in the thymus [46]. In addition, it drives epigenetic remodelling of CD8 + T cells towards an exhausted phenotype in chronic infections [47,48]. Notably, in the CD4 + T cells, we identified a hypermethylation in cases at the intergenic region near the TSS of the TOX gene. The methylation at this region had a weak negative correlation (−0.157) with gene expression, suggesting that hypermethylation of this locus in cases may favour T cell differentiation towards an effector phenotype in type 1 diabetes patients. It was reported that autoreactive CD8 + T cells with an exhaustion-like profile associated with slow type 1 diabetes progression [49] and preservation of beta cell function in type 1 diabetes patients [50]. Further, genetic variants at TOX are among novel risk loci for type 1 diabetes [51].
Our recent study aimed at identifying umbilical cord blood DNA methylation associated with type 1 diabetes progression [19], as well as a previous study by another group [12], did not find any difference in DNA methylation, suggesting that significant DNA methylation differences in children progressing to the disease are not yet present at birth or that they are too subtle to be identified with the current methods and cohort size.
A recent study by Johnson et al examined DNA methylation changes in blood samples longitudinally collected before diagnosis of type 1 diabetes from children enrolled in the prospective DAISY cohort using Illumina 450 K and EPIC platforms [14]. They observed 28 DMRs associated with later development of the disease. We found no overlap with the CpG sites identified in the study by Johnson et al. This may perhaps be due to different methodologies and to the fact that Johnson et al analysed peripheral whole blood samples, whereas we examined DNA methylation in purified cell fractions of PBMCs. Despite these differences, both studies identified pre-seroconversion changes in DNA methylation, suggesting that these may predict disease progression before the appearance of autoantibodies. However, more studies are needed to validate the findings in independent, larger cohorts, with replication in other populations.
Interestingly, Paul et al [12] identified one DMC (hg19: chr10:74058002) in purified CD4 + T cells between individuals with type 1 diabetes and their healthy co-twins in the intergenic region between DDIT4 and DNAJB12. In our study, we found a DMR (chr10:74082048--74,082,090) region hypermethylated in cases in CD4 + T cells in preseroconverted samples located in the same intergenic region (Fig. 6c).
In conclusion, our results provide novel insights into cell type-specific DNA methylation changes associated with type 1 diabetes development. These findings provide the basis for further studies to find an early methylation signature, which may be useful for disease prediction and management.
Genomics Centre (FFGC), supported by Biocenter Finland, for their assistance. We also acknowledge the Finnish Centre for Scientific Computing (CSC), Elixir Finland and the Aalto Science-IT project for computational resources. The graphical abstract was created with BioRender.
Data availability The analysed datasets are available through ArrayExpress accession E-MTAB-11088. Due to privacy reasons, the raw data are available through the corresponding author upon reasonable request.
Funding Open Access funding provided by University of Turku (UTU) including Turku University Central Hospital. R Lahesmaa was supported by the Academy of Finland (AoF) (grants 292335, 294337, 319280, 31444, 319280, 329277, 331790) and Business Finland, and grants from JDRF, the Novo Nordisk Foundation (grant NNF19OC0057218), the Sigrid Jusélius Foundation (SJF), the Jane and Aatos Erkko Foundation, the Finnish Diabetes Foundation and the Finnish Cancer Author's relationships and activities The authors declare that there are no relationships or activities that might bias, or be perceived to bias, their work.
Contribution statement IS and UUK contributed to data analysis, interpreted the results and wrote the manuscript. EL was responsible for data analysis and contributed to writing. TG and RK were responsible for pyrosequencing analysis. TH and VT played a key role in the DIABIMMUNE study and provided clinical data and samples. SJ, AL and VS participated in data analysis and interpretation. HK contributed to the sample and data acquisition. R. Lund, LLE and HL contributed to data analysis and interpretation. R Lund also participated in supervision of the laboratory work, provided expertise in the bisulphite sequencing. LLE supervised SJ, AL and VS for data analysis and participated in the interpretation of the results. HL supervised EL for data analysis and participated in the interpretation of the results. MK led the DIABIMMUNE study, provided samples and clinical data, initiated the study with R. Lahesmaa, interpreted the results and edited the manuscript. R. Lahesmaa initiated and designed the study setup, supervised the study, interpreted the results and wrote the manuscript. All authors approved the final version. R. Lahesmaa and EL are the guarantors of the study, accept full responsibility for the work, had access to the data and controlled the decision to publish.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.