Longitudinal changes in DNA methylation during the onset of islet autoimmunity differentiate between reversion versus progression of islet autoimmunity

Background: Type 1 diabetes (T1D) is preceded by a heterogenous pre-clinical phase, islet autoimmunity (IA). We aimed to identify pre vs. post-IA seroconversion (SV) changes in DNAm that differed across three IA progression phenotypes, those who lose autoantibodies (reverters), progress to clinical T1D (progressors), or maintain autoantibody levels (maintainers).

Methods: This epigenome-wide association study (EWAS) included longitudinal DNAm measurements in blood (Illumina 450K and EPIC) from participants in Diabetes Autoimmunity Study in the Young (DAISY) who developed IA, one or more islet autoantibodies on at least two consecutive visits. We compared reverters - individuals who sero-reverted, negative for all autoantibodies on at least two consecutive visits and did not develop T1D (n=41); maintainers - continued to test positive for autoantibodies but did not develop T1D (n=60); progressors - developed clinical T1D (n=42). DNAm data were measured before (pre-SV visit) and after IA (post-SV visit). Linear mixed models were used to test for differences in pre- vs post-SV changes in DNAm across the three groups. Linear mixed models were also used to test for group differences in average DNAm. Cell proportions, age, and sex were adjusted for in all models. Median follow-up across all participants was 15.5 yrs. (interquartile range (IQR): 10.8-18.7).

Results: The median age at the pre-SV visit was 2.2 yrs. (IQR: 0.8-5.3) in progressors, compared to 6.0 yrs. (IQR: 1.3-8.4) in reverters, and 5.7 yrs. (IQR: 1.4-9.7) in maintainers. Median time between the visits was similar in reverters 1.4 yrs. (IQR: 1-1.9), maintainers 1.3 yrs. (IQR: 1.0-2.0), and progressors 1.8 yrs. (IQR: 1.0-2.0). Changes in DNAm, pre- vs post-SV, differed across the groups at one site (cg16066195) and 11 regions. Average DNAm (mean of pre- and post-SV) differed across 22 regions.

Conclusion: Differentially changing DNAm regions were located in genomic areas related to beta cell function, immune cell differentiation, and immune cell function.

1 Introduction

T1D is an autoimmune disorder with significant long-term morbidity. The pre-clinical phase is defined by the appearance of autoantibodies against pancreas cell antigens, termed islet autoimmunity (IA). There is strong evidence to support autoantibodies as a biomarker of T1D risk (1). However, IA is dynamic. While progression to T1D or multiple autoantibodies has been well characterized, a subset of individuals lose autoantibody positivity (2) and revert back to an autoantibody negative state. Autoantibody reversion was first described by Spencer et al (3) in a cohort of 685 individuals with a first degree relative affected by T1D. After 5 years, 7/20 developed T1D, 1 remained AB positive and 12/20 reverted. Transient autoantibody positivity has been described in several additional studies (4–6). However, these historical studies describing the transient nature of autoantibodies are difficult to interpret due to the development of more accurate autoantibody tests as well as differences in the definition of reversion. Vehik et al (2) conducted the most comprehensive and rigorous study of reversion in current literature. Among 596 individuals enrolled in The Environmental Determinants of Diabetes in the Young (TEDDY) study who developed one or more persistent autoantibodies, 21% reverted to an antibody negative state. Seroreversion was associated with significantly decreased risk of T1D (hazard ratio: 0.14, 95% CI: 0.04-0.59). Understanding the unique protective mechanisms occurring prior to or following IA that lead to IA reversion may have important implications for development of interventions that delay or prevent progression to T1D.

Genetic variation is a well-established risk factor for T1D (7). However, heterogeneity in disease concordance among monozygotic twins (8) as well as temporal changes in both T1D incidence (9) and age at T1D onset (10) in population studies have created a strong interest in the role of the environment in the etiology of T1D. Epigenetic modifications such as DNA methylation (DNAm) may represent a mechanistic pathway between genetic susceptibility, environmental exposures, and progression or reversion of IA. Epigenetics broadly describes a class of modifiable mechanisms that can regulate gene expression and are sensitive to external stimuli (11). DNAm is a frequently studied epigenetic biomarker that is postulated to play a role in autoimmune diseases as epigenetic mechanisms are important regulators of immune cell differentiation, plasticity and function (12, 13). DNAm changes prior to and during the IA phase may provide key information about underlying epigenetic profiles that explain progression or reversion from IA.

Previous epigenome wide studies have identified significant associations between DNAm and T1D (14–17). However, associations have been inconsistent and many of the studies have focused on static and/or post-T1D differences in DNAm between cases and controls (14–16). Although important in understanding the etiology of T1D, DNAm differences obtained from a single time point are difficult to interpret as it is not possible to determine when the changes occurred and moreover, whether they are the cause or consequence of the disease process. Understanding the timing of the changes is key to identifying external factors that cause these changes and therefore, may be amenable to preventative interventions. The purpose of this study was to test DNAm obtained before and after IA seroconversion (SV) in the Diabetes Autoimmunity Study in the Young (DAISY). We aimed to identify pre vs. post-SV changes in DNAm that differed across three distinct IA progression phenotypes, those who lose autoantibodies (reverters), progress to clinical T1D (progressors), or maintain autoantibody levels (maintainers).

2 Materials and methods

2.1 Study population

We reviewed individuals from the Diabetes Autoimmunity Study in the Young (DAISY) who developed islet autoimmunity (IA) between February 1994 and February 2019. DAISY is a longitudinal birth cohort study that includes n=2544 children at high risk for T1D. Subjects are recruited from two high risk populations, those with a first degree relative (FDR) with T1D or those with a high-risk genotype, [defined as DRB1*04, DQB1*0302/DRB1*0301, DQB1*0201 (DR3/4 DQ8)]. Subjects complete study visits at 9, 15, and 24 months. Following the 24-month visit, study visits occur annually. As described previously (18), radio-immunoassays were used to test serum samples for autoantibodies to insulin (IAA), GAD65 (GAA), and IA-2 (IA-2A). Prior to 2010, GADA and IA-2A were tested using a combined radioassay (19). The National Institute of Diabetes and Digestive and Kidney Diseases harmonized assay was used to test for GADA and IA-2A after 2010 (20). Serum samples from individuals positive for GAD65, IAA, or IA-2 were tested for ZnT8A following development and implementation of the ZnT8 assay (21). If autoantibodies are detected, participants return for study visits every 3-6 months.

Islet autoimmunity (IA) was defined as the presence of one or more autoantibodies (see above) on at least two consecutive visits 3-6 months apart. The first visit among these consecutive autoantibody positive visits designated the start of IA, referred to as seroconversion (SV) throughout the remainder of the manuscript. We defined the three autoimmune progression phenotypes based on the autoantibody testing. The reverter group was defined as individuals who reverted for all autoantibodies during two or more consecutive visits, did not develop T1D, and were autoantibody negative for all autoantibodies at their last DAISY visit. The maintainer group was defined as individuals who continued to test positive for islet autoantibodies and did not develop T1D at the time of their last visit. The progressor group was defined as individuals who developed clinical T1D.

Among individuals who developed IA during DAISY and underwent autoantibody testing for a minimum of two or more study visits (n=213), we excluded individuals for the following: missing a pre- or post-SV blood sample (n=54), onset of IA unclear due to gaps (>365 days) in study visits (n=2), missing study visit prior to initial pre-SV positive visit (n=14). The Colorado Multiple Institutional Review Board approved all DAISY protocols (COMIRB 92-080). Informed consent and assent, if appropriate, was obtained from the parents/legal guardians of all children prior to participation in any research related activities.

2.2 Methylation measurements

Methylation measurements were obtained from peripheral whole blood samples collected at multiple time-points in individuals from DAISY. The Infinium HumanMethylation 450K Beadchip platform (Illumina, San Diego, CA, USA) was used to obtain methylation measurements on a subset of samples. The 850K Infininium MethylationEPIC BeadChip (Illumina, San Diego, CA, USA) was used to obtain measurements on the remaining samples. Two platforms were used due to changes in technology during the course of the study. Samples were randomly assigned to the two platforms making sure all timepoints from the same individual were included on the same platform.

DNA was bisulfite converted using the Zymo EZ DNA Methylation kit (Zymo Research, CA, USA). The bisulfite-converted DNA was labeled with fluorescent dyes and hybridized to 450K and 850K DNAm arrays. Samples were arranged on the plates in a specific sequence to minimize within and between batch effects (plate effects are represented by first 11 digits of the array variable on GEO). The minfi (v1.12.0) package (22) in R (v3.5.2) was used to perform quality control (QC) checks at the sample level. The processing pipeline is described in greater detail in Vanderlinden et al (23).

The DNAm probes were annotated to the genome based on the hg19 genome build using the Illumina annotation manifest files. Non-autosomal CpGs or CpGs located within or near (<2 base pairs) known single nucleotide polymorphisms (SNPs) were excluded. CpG sites with a beta range <3% on both platforms were removed from analysis. A total of n=198,008 overlapping DNAm probes met our filtering criteria and were used in subsequent analyses. Normalized M-values (SeSAMe (v1.0.0) pipeline with Noob normalization) were used in all statistical analyses. We use the term DNAm probe and the probe identifier when referring to the data in the Methods and Results. However, each probe is designed to measure DNAm at a single CpG site which is used as a more general term in the Discussion. See Figure 1 for an overview of the study methods.

Figure 1

Figure 1 Summary of methods used to identify and prioritize DNAm candidates. Description: We used an epigenome wide association study design to identify differentially methylated positions (DMP) associated with the three islet autoimmunity progression phenotypes, reverters, maintainer, or progressors. We used two DMP models (1) an interaction model that tested whether changes in DNA methylation (DNAm) levels at single CpGs pre-IA versus post-IA differed across groups and (2) a group effect model that tested whether average methylation levels (pre- and post-IA) differed across groups. We also performed regional analyses (differentially methylated regions or DMRs) based on single CpG sites from the two models to identify regions with consistent methylation effects. We identified regions where average regional methylation levels differed between groups (μDMRs) as well as regions where changes in regional methylation levels pre- vs post-IA differed across groups (ΔDMRs). In order to prioritize regions, we tested whether the DNAm candidates identified in our analysis were associated with gene expression levels post-SV, an expression quantitative trait methylation analysis (eQTM). To account for the multiple CpGs within each DMR, we used a principal component analysis to capture common patterns across all CpGs included in the DMR. We identified cis-eQTMs (midpoint of region +/- 500 KB of the TSS of the gene) by testing the correlation between gene expression and the 1st principal component. We also tested the correlation between DNAm candidates and metabolite levels obtained from overlapping samples, a metabolite quantitative trait methylation analysis (metQTM). We used a principal component analysis to capture common patterns across all CpGs included in the candidate DMRs. We tested the correlation between metabolite levels and the 1st principal component. CpGs are represented by lollipop plots in the figure.

2.3 Overlapping gene expression measurements

Gene expression data were available in a subset of individuals (n=36) at the post-SV visit. RNA processing and quantification is described in greater detail in Carry et al (24). In brief, paired end sequencing was performed using the Illumina NovaSEQ 6000™ system and samples were quantified against the Ensembl reference transcriptome (hg19, version 87) using the RSEM algorithm (25). Data were quantile normalized using DESeq2 (26), re-normalized using RUV (27), and then transformed using the regularized log function (26). The transformed data were used in all subsequent statistical analyses.

2.4 Overlapping metabolomics measurements

Untargeted metabolomics data were available in a subset (n=110) of individuals at both the pre-SV and post-SV visits. Metabolomics processing and quantification is described in greater detail in Carry et al (28). In brief, non-fasting plasma samples were used to quantify metabolite levels using three untargeted panels, HILIC panel: HILIC-QTOF MS/MS (29), GCTOF panel: GC-TOF-MS (30), and Lipid panel: CSH-QTOF MS/MS (31). BinBase (32) was used to process and annotate the GC-TOF-MS data. MS-Dial (33) was used to process and annotate the liquid chromatography (LC), CSH-QTOF-MS and HILIC-QTOF-MS, data. LipidBlast (34) and Massbank of North America were also used to annotate the complex lipids (http://mona.fiehnlab.ucdavis.edu/). Metabolomic data were normalized using the systematic error removal using random forest (SERRF) algorithm (35). All metabolites were Box-Cox transformed prior to statistical analysis.

2.5 Genetic ancestry

Ancestry principal components (PC) were estimated for all study participants from genetic data collected in DAISY. Sample processing and genotyping were performed at the University of Virginia School of Medicine Center for Public Health Genomics based on exome genotyping (Illumina HumanCoreExome-24 BeadChip, N=283) or whole genome sequencing (N=162) from the larger DAISY population, see Buckner et al (36) for a more complete description of the genetic processing and calculation of the genetic ancestry PCs.

2.6 Statistical analyses

The overall methods workflow is summarized in Figure 1. Linear mixed models were used to test for differences in DNAm between the pre- and post-SV visit across reverters, maintainers, and progressors (autoimmune phenotype*visit interaction). Separate linear mixed models were also used to test for differences in average DNAm (mean of the DNAm levels at the pre- and post- SV visits) between the autoimmune phenotypes (group effect). Platform (EPIC vs 450K), age, sex, and cell proportions (estimated using the minfi (v1.12.0) package (22) implementation of the Houseman method) were adjusted for in all models. The group effect models were also adjusted for population ancestry (see Supplementary Material for complete description of ancestry data). Ancestry data (1^st 2 PCs) were unavailable for 2 individuals in the group effect model and thus, these individuals were not included in this analysis. See Appendix 1 (Data Sheet 1) for the linear mixed model code. We did not adjust for ancestry in the interaction (autoimmune phenotype*visit) models because the interaction models test for within individual differences, and thus are less likely to be impacted by time invariant confounders such as population ancestry. The Benjamini Hochberg false discovery rate (FDR), was used to correct for multiple comparisons (37). Significance was assessed based on the FDR adjusted p-value <0.10. Model diagnostics are described in the Supplementary Files (Data Sheet 2), see Appendix 2, Figures A–C and Table A.

The comb-p python software package (38) was used to identify differentially methylated regions (DMRs). Within the comb-p pipeline, we used a seed p-value of 0.1 and then searched for adjacent probes within a window of 500 bases, using a step size of 50 bases. Comb-p combines probes within this window and then calculates an overall, spatially corrected p value for the entire region based on the Stouffer-Liptak method. The Sidak method is used to adjust the overall regional p values for multiple testing. Regional analyses were performed based on the individual DNAm probes from the interaction (post- vs pre-SV changes by autoimmune phenotype), referred to as differentially changing DMRs (ΔDMR) throughout the remainder of the manuscript. Regional analyses were also performed based on DNAm probes from the main effect model (differences in average of pre- and post-SV DNAm between groups), referred to as average DMRs (μDMR) throughout the remainder of the manuscript. For both regional analyses, we reviewed all regions with ≥4 DNAm probes that were significant at the combined Sidak adjusted region p value of 0.10. Because the interaction and group effect p values are based on a two degree of freedom test (numerator degrees of freedom for the overall F-test), it is possible for the DMR to capture a set of DNAm probes with similar p values but substantial heterogeneity in the directions of effect within the three groups. Therefore, for the ΔDMRs, we retained regions with a consistent direction of effect, defined as a region where the direction of change in DNAm between the two visits (hyper methylation or hypo methylation) was consistent across 100% of the DNAm probes within the region in one or more of the study groups. For the μDMRs, we retained regions where the direction of effect (hypo or hypermethylation) for one or more of the pairwise group comparisons was consistent across 100% of the DNAm probes included in the region.

2.7 Expression quantitative trait methylation analysis: correlation between gene expression and DNAm candidates

In order to better understand our primary DNAm results, we tested the correlation between gene expression levels and our DNAm candidates, one DMP, 11 ΔDMRs, and 22 μDMRs in a subset of individuals (n=36, see Appendix 3, Table B) with methylation data pre- and post-SV as well as gene expression data post-SV. First, linear mixed models were used to regress out age, sex, platform, and cell proportions from the DNAm values at each of the candidate CpG sites. Ancestry PC1 and ancestry PC2 were also regressed out from all CpG sites included in the μDMRs candidate regions. Next, using the residuals from the linear mixed models, the within individual differences in DNAm (post-SV minus pre-SV) were used to represent changes in DNAm between the study visits for each of the CpG sites included in the ΔDMRs. The average residual values from the post-SV and pre-SV study visits were used to represent average methylation for each of the CpG sites within the μDMRs. Next, we performed a principal component analysis of DNAm levels across the region-specific CpG sets. For each DMR, the first PC was extracted for subsequent testing, allowing us to consider all CpG sites together rather than testing many individual sites separately. Linear regression models were then used to regress out the effects of age and sex from the gene expression levels. Finally, Spearman correlation coefficients were used to test the correlation between DNAm and gene expression residuals. We looked for cis-eQTMs, defined as genes significant at the FDR adjusted p value of 0.10 where transcription start site was +/- 500 KB of the midpoint of the DMR. FDR adjustment was based on the total number of DNAm cis-gene pairs (256 transcript DNAm pairs for the ΔDMR candidates and 544 transcript DNAm pairs for the μDMR candidates).

2.8 Metabolite quantitative trait methylation analysis: correlation between metabolite levels and DNAm candidates

We tested the correlation between DNAm and untargeted metabolite levels in a subset of our study population (n=110, see Appendix 3, Table B) with DNAm and metabolomics data available both pre- and post-SV. Only data from overlapping samples was included in this supplementary analysis. Linear models were used to regress age and sex from the Box-Cox transformed metabolite levels at each visit. Consistent with the DNAm methods, using the residuals from the linear mixed models, the difference between metabolite residuals at each visit (post-SV minus pre-SV residuals) was used to represent change in metabolites and the average residual values (average of post-SV and pre-SV residuals) were used to represent average metabolite values. For the ΔDMR candidates and the single DMP candidate, linear regression models were then used to test the correlation between the change in metabolites versus the ΔDMR PCs (described above) as well as the single DMP candidate. For the μDMR candidates, linear regression models were then used to test the correlation between average metabolite levels versus the μDMR PCs (described above). False discovery (FDR) rate adjusted p values were calculated for all individual metabolite DNAm candidate pairs according to methods described by Benjamini and Hochberg (37). FDR adjusted p values were calculated separately for each platform. Only annotated metabolites from the HILIC (81 metabolites), Lipid (373 metabolites), and GC-TOF (98 metabolites) panels were evaluated in subsequent analyses. Metabolites were evaluated at an FDR adjusted p value of 0.10.

3 Results

3.1 Study population

The final study population included 60 individuals in the maintainer group, 42 individuals in the progressor group, and 41 individuals in the reverter group. At both the pre-SV and post-SV visits, age differed by group, and the estimated cell proportions differed by group at the post-SV visit (Table 1). At the time of data analysis, duration of follow-up, defined as median time from the initial visit to the development of T1D or last study visit, was 9.3 years (IQR: 6.1 to 12.3 years) for the progressors, 16.5 years for the maintainers (IQR: 14.3 to 20.9 years) and 16.6 years for the reverters (IQR: 15.2 to 20.2 years).

Table 1

Table 1 Demographics and clinical characteristics.

The specific autoantibody subgroups present at the onset of seroconversion in the three groups are described in greater detail in Appendix 4 (Data Sheet 4), Table C. As expected, the prevalence of multiple autoantibodies at serconversion was higher in progressors (31%) relative to maintainers (18%) and reverters (0%). Across the entire islet autoimmunity follow-up period, the occurrence of multiple autoantibodies at one or more study visit(s) following IA seroconversion was also higher in progressors (86%) compared to maintainers (58%). Among reverters, 10% developed multiple autoantibodies at one of more study visit(s) during the time period between seroconversion (IA onset) and seroreversion.

3.2 Differentially methylated position analysis

Change in methylation at the DNAm site cg16066195 on chr 7 was significantly (FDR adjusted p value=0.0174) different across groups. The reverter group was characterized by an increase in DNAm between pre- and post-SV visits (ie, a positive slope) whereas the progressor and maintainer groups were characterized by no change or a decrease in DNAm (Figure 2). This site is an island CpG site (CpG island chr7:73703458-73704127) that maps to an area near the CLIP2 gene.

Figure 2

Figure 2 Changes in DNAm between the pre- and post-SV visits at cg16066195 across the three IA progression phenotypes. Description: (A) provides the average methylation M-values and corresponding 95% confidence intervals within the three IA progression phenotypes pre- and post-SV. (B) describes the individual level changes in methylation m-values (y-axis) between the post-SV visit relative to the pre-SV visit in the three IA progression phenotypes (x-axis). Positive values represent increasing DNAm whereas negative values represent decreasing methylation between visits. All DNAm values in (A, B) have been adjusted for age, sex, and cell proportions.

We also tested whether average DNAm (mean of DNAm levels pre- and post-SV) differed across groups. No DNAm probe was significant at the FDR adjusted alpha level of 0.10.

3.3 Differentially methylated region analysis

We also tested for genomic regions (Figure 1). In contrast to the single CpG site (DMP) analysis, the regional analysis allowed us to identify multiple CpG sites that demonstrated similar DNAm changes between the pre- and post-SV visits across the three study groups (ΔDMRs). We focused on FDR significant regions of ≥4 DNAm probes where the direction of the change in DNAm (between the pre-SV and post-SV visits) was consistent (100% of probes changed in a similar direction) within one or more of the groups. We identified 11 candidate DMRs (Table 2; Figure 3).

Table 2

Table 2 Regions where DNAm changes between the post- and pre-SV visits were consistently different across groups (group*visit interaction).

Figure 3

Figure 3 Differentially changing methylation region on chromosome 20 where changes in DNAm (pre- vs post-SV) differed across the three IA progression phenotypes. Description. Region on chromosome 20 loc 57426538 to 57427974 (ΔDMR1) where the change in DNA methylation (DNAm) post- vs pre-SV differed across groups. In the top panel, each dot represents the within group slopes (y-axis) or changes in DNAm m-values between the post-SV and pre-SV visit at each of the CpG sites included ΔDMR 1. The x-axis represents the position (mb) of the CpGs within the region. All slope values were adjusted for age, sex, and cell proportions. Positive values indicate methylation values increased following IA seroconversion whereas negative values indicate methylation decreased following IA seroconversion. The dashed lines represent the average slope value within each group across the entire region. The middle panel represents the location of the region (black solid square) relative to the closest genes, GNAS and ATP5E (red solid boxes). There are multiple known isoforms for GNAS and ATP5E, the bottom panel displays the most biologically relevant or consensus transcript based on the Ensembl database. The red line on the ideogram, bottom of the figure, represents the location of GNAS and ATP5E on chromosome 20.

We also tested for regions where the average DNAm levels at the pre- and post-SV visits differed across the groups (μDMRs). We identified 22 FDR significant μDMRs of ≥4 DNAm probes where the direction of the pairwise group differences in DNAm was consistent across all CpG sites included in the region (Table 3; Figure 4).

Table 3

Table 3 Regions where average of post- and pre-SV DNAm levels were consistently different across groups (group main effect).

Figure 4

Figure 4 Differentially methylated region on chromosome 12 where average (pre- and post-SV) methylation levels differed across the three IA progression phenotypes. Description. Region on chromosome 12 loc 2943902 to 2944481 (μDMR4) where average DNA methylation (DNAm) levels, post- and pre-SV, differed across groups. In the top panel, each dot represents the average DNAm value (y-axis) at each of the CpG sites included μDMR4. The x-axis represents the position (mb) of the CpGs within the region. All DNAm values were adjusted for age, sex, cell proportions, and genetic ancestry. The dashed lines represent the average methylation value within each group across the entire region. The middle panel represents the location of the region (black solid square) relative to the closest genes, ITFG2 and NRIP2 (red solid squares). There are multiple known isoforms for ITFG2 and NRIP2, the figure displays the most biologically relevant or consensus transcript based on the Ensembl database. The red line on the ideogram, bottom of the figure, represents the location of ITFG2 and NRIP2 on chromosome 12.

3.4 eQTM candidate prioritization

We tested the correlation between DNAm and cis- gene expression levels in a subset of overlapping samples. The availability of individual level DNAm data allowed us to look at the entire DMR together. Based on the ΔDMR candidates, we identified two FDR significant cis eQTMs representing one DMR and two gene transcripts, GNAS and ATP5E (ΔDMR1, region on chromosome 20, see Table 4). Within this region, increased DNAm post- vs pre-SV was positively associated with expression of GNAS and ATP5E (see Table 4).

Table 4

Table 4 Summary of FDR significant cis-eQTMs representing correlation between differentially changing methylation regions and gene expression post- SV.

3.5 Metabolite quantitative trait methylation analysis candidate prioritization in overlapping samples

We tested whether the single DMP candidate, cg16066195, as well as the candidate DNAm regions identified in our primary analysis were associated with metabolite levels. Consistent with the eQTM analysis, we regressed out age and sex from annotated metabolites and then tested the correlation between annotated metabolites versus DNAm regional PCs. Based on the ΔDMR candidates, we identified 26 annotated metabolites from the Lipid panel that were correlated with 4 DMRs (see Table 5; Figure 5). ΔDMR 8 was correlated with multiple lipids, primarily PCs, ΔDMR 5 was also correlated with multiple lipids, primarily correlated with TGs (fats). ΔDMR 9 and ΔDMR 2 were correlated with a single lipid, an ether lipid, and a TG, respectively. Metabolite candidates primarily consisted of odd-chain fatty acid containing lipid species (OCFA). Furthermore, the majority of the metabolites (29/30) were positively correlated with increasing DNAm levels. The μDMR candidate regions as well as the single DMP candidate were not significantly associated with metabolite levels at our FDR adjusted cutoff of 0.10.

Table 5

Table 5 Secondary metQTM analysis of the association between pre- versus post-SV change in methylation across the DMRs and pre- versus post-SV change in metabolite levels.

Figure 5

Figure 5 Differentially changing region on chromosome 6 (post- vs pre-SV) that was positively correlated with changes in lipid metabolites (post- vs pre-SV). Description: Region on chromosome 6 loc 170597377 to 170597899 (ΔDMR8) where the change in DNA methylation (DNAm) post- vs pre-SV differed across groups. In the top left (A), each dot represents the within group slopes (y-axis) or changes in methylation m-values between the post-SV and pre-SV visit at each of the CpG sites included ΔDMR 8. The x-axis represents the position (mb) of the CpGs within the region. Positive values indicate DNAm values increased following IA seroconversion whereas negative values indicate DNAm decreased following IA seroconversion. The dashed lines represent the average slope value within each group across the entire region. The top right (B) represents the association between DMR wide DNAm captured by the 1^st PC (x-axis) and changes in metabolite values (y-axis) between the post- and pre-SV visits. DNAm and metabolite expression values have been standardized to facilitate the interpretation of the slope as a 1 standard deviation increase in the change in metabolite levels between the post- and pre-SV visits per 1 standard deviation increase in the change in methylation between post- and pre-SV visits. The bottom panels (C, D) represent the average metabolite levels and corresponding 95% confidence intervals within the three groups pre- and post-SV. All DNAm and metabolite values were adjusted for age, sex, and cell proportions.

4 Discussion

Epigenetic biomarkers are appealing in the study of complex diseases such as T1D based on their heritability, role in gene expression, and responsiveness to external stimuli. Epigenetic effects in observational studies are challenging to interpret because it is often not possible to determine whether DNA methylation (DNAm) is causative or secondary to the disease process. A strength of our study is the longitudinal analysis of DNAm levels both before and after the onset of IA. We identified a single CpG site as well as genomic regions where changes in DNAm between the post-SV and pre-SV visits were significantly different across the IA progression phenotypes. We also identified regions where average DNAm levels pre- and post-SV differed across the progression phenotypes. Together, the DNAm regions have potential biological relevance to T1D etiology based on their potential role in immune and beta cell function.

We identified a DNAm site, cg16066195, on chromosome 7 where DNAm levels increased between the pre- and post-SV visits among individuals who reverted to an IA negative state (reverters) compared to progressors (who showed no change in DNAm) and maintainers (who showed decreasing DNAm, Figure 2). This island CpG is located near the transcription start site for the protein coding gene CLIP2. In a mouse model of diet induced changes in beta cell expression, CLIP2 gene expression was significantly downregulated among mice fed a carbohydrate containing diabetogenic high-fat diet relative to mice fed a diabetes-protective carbohydrate free high-fat diet (39). Furthermore, SNPs within CLIP2 (rs2528994 and rs512023) have demonstrated modest associations with T2D in both the Diabetes Genetics Initiative (40) and the Wellcome Trust Case Control Consortium (41).

Our methylation analysis also identified numerous regions where average methylation post- and pre-SV differed across the autoimmune phenotypes in areas of the genome potentially relevant to T1D etiology. We identified a DMR on chromosome 12, μDMR4, characterized by hypermethylation in the reverter group relative to the progressor and maintainer groups (Figure 4). This includes 4 probes that, based on the ENCODE Project Consortium (42), are located in a known enhancer region. Three of the four probes within this region are located within the transcription start site for NRIP2, predicted to act upstream or within the notch signaling pathway (43). This pathway is relevant to T1D (44) based on its role in immune cell differentiation and function (45) as well as pancreas development (46), islet cell function (47), and islet cell survival (48). All four probes within μDMR4 are also located within the 5’UTR region for ITFG2, a gene expressed in numerous tissues including immune cells. Mouse and in vitro models have demonstrated that ITFG2 deficiency alters B cell maturation and migration (49). In a lupus mouse model, MRL/lpr, autoimmunity development occurred earlier and was more severe in ITFG2 deficient mice (49). Together, these findings suggest a potential role for ITFG2 in B cell differentiation and as a potential regulator of autoimmunity. Although, average methylation within DMR4 was not correlated with expression of ITFG2 or NRIP2 in our secondary eQTM analysis, three probes within μDMR4 (cg05194726; cg06997549; cg02852959) were correlated with expression of both ITFG2 and NRIP2 in whole blood based on the BIOS QTL browser (50), an online resource that provides a searchable database of FDR significant associations between DNAm and gene expression (eQTM). Additional work is needed to understand the connections between methylation within this region on chr 12, ITFG2 expression, NRIP2 expression, and T1D etiology.

We also identified several regions of differentially changing DNAm that are potentially relevant to T1D etiology based on known associations between DNAm in these regions and relevant environmental risk factors. We identified a region on chr 20 near the GNAS/GNASAS loci, ΔDMR 1, that was characterized by decreasing DNAm pre- vs post-SV in maintainers and progressors relative to reverters (Table 2; Figure 3). Based on the ENCODE Project Consortium (42), 25 of the 29 probes in ΔDMR 1 are located within a DNAase hypersensitivity region and 4 probes are known to interact with transcription factor binding. DNAm in this region is responsive to environmental stressors. Umbilical cord blood DNAm near GNAS was altered among infants born to a mother affected by gestational diabetes (GDM), a disorder characterized by glucose intolerance during pregnancy (51). Based on the Dutch Hunger Winter Families Study (52), siblings exposed to the war-time Dutch Hunger Winter famine were associated with persistent changes in DNAm in a region near the GNASAS locus relative to their unexposed siblings (53). The direction and magnitude of effect depended on timing of exposure and sex of the exposed individual (53). DNAm among exposed siblings was also altered near another gene implicated in metabolic disease MEG3 (53), a gene that mapped to ΔDMR4 which was also characterized by decreasing methylation among progressors and maintainers relative to reverters (Table 2). Interestingly, both the GNAS (54) and MEG3 (55) genes are maternally imprinted. Loss of maternal imprinting should be investigated as a potential mechanism in the etiology of T1D using whole-genome bisulfite sequencing in order to provide a higher density representation of DNAm changes within imprinted areas of the genome.

The secondary eQTM analysis in a subset of overlapping samples confirmed that changes in methylation within ΔDMR1 were associated with expression of GNAS. Increased methylation post- versus pre-SV was associated with higher levels of GNAS expression at the post-SV visit in a subset of overlapping samples. GNAS is an important regulator of insulin secretion in beta cells (56). GNAS silencing results in decreased insulin secretion and insulin content (56). GNAS encodes the G protein subunit alpha which also plays a role in the interaction between antigen presenting cells and T helper cell differentiation (57). Mice with dendritic cells deficient for GNAS result in a phenotype characterized by preferential Th2 differentiation, Th2 type inflammation, and subsequent development of allergic asthma (57). Overlap between autoimmunity and atopic conditions have long been hypothesized based on disruptions in similar immune pathways (58). Positive associations between childhood asthma and subsequent T1D development have been observed in several countries (59–61). Overall, our results suggest that maintenance of DNAm levels near GNAS during IA may represent a unique protective mechanism in reverters.

In order to further characterize the DNAm regions identified in the primary analysis, we tested the correlation between changes in DNAm and changes in annotated metabolites (metQTM). Four differentially changing DMRs were correlated with changes in 26 unique lipid metabolites (Table 5). ΔDMR 8, characterized by increasing methylation in progressors (Figure 5), was correlated with 18 of the 26 lipid metabolites. This region of differentially changing methylation is notable based on its location in an open chromatin region within the body of the DLL1 gene on chr. 6. As a notch signaling ligand, DLL1 controls the differentiation of pancreatic progenitor cells into exocrine versus endocrine cells (46). The loss of DLL1 results in early progenitor cell differentiation and an overabundance of endocrine cells (46). A recent mouse model confirmed DLL1 is also relevant to islet cell function in the mature pancreas based on its high level of expression in beta cells and corresponding role in insulin secretion (47). Furthermore, DLL1 plays an important role in differentiation of B cells and the development of antigen secreting cells; the presence of DLL1 influences AB titer levels and isotype switching (45). Additional work is needed to understand the connection between a concordant increase in lipid levels and DNAm within the DLL1 gene following seroconversion.

Our secondary metQTM was unique in that DNAm and metabolite levels were available pre- and post-SV in a subset of overlapping samples. This analysis revealed a consistent positive association between increasing lipid metabolite levels, post- vs pre-SV, and increasing DNAm levels across several regions (25 of the 26 unique lipid metabolites were positively correlated with DNAm changes, see Table 5). Numerous studies (62–68) have reported associations between dysregulation in lipid levels and T1D. Although lipid levels have been shown to be influenced by age at sample collection/timing of sample collection relative to onset of IA and type of first appearing autoantibody, prior research suggests lower lipid levels, including sphingomyelins and phosphatidylcholines, are generally associated with increased risk of T1D and/or IA (62–68). In our study, increasing lipid levels, in particular phosphocholines, following the onset of IA were strongly correlated with increasing methylation within ΔDMR8. This region was characterized by increasing methylation within the progressor group. However, as demonstrated in Figure 5, the lipid metabolite most strongly correlated with DNAm changes in this region, Phosphatidylcholine (35:4), was lower in the progressor group prior to SV and then subsequently increased following the onset of IA, suggesting higher levels of lipids within the progressor group may be unique to changes that occur following seroconversion.

There was a high prevalence of odd-chain fatty acid (OCFA) containing lipid species among the metabolites correlated with DNAm changes. Recently, there has been increased recognition of OCFA in plasma and their potential biological relevance (69). OCFA levels have been associated with glucose homeostasis, insulin resistance, T2D, and BMI (69, 70). Pfleuger et al (71) observed higher levels of odd-chain triglycerides among autoantibody positive versus negative children in BABYDIAB. This parallels the concordant post-seroconversion increase in OCFA levels and DNAm near the DLL1 gene (ΔDMR 8) among progressors (Figure 5) in the current study. OCFA have been proposed a marker of dairy intake which has been positively correlated with progression to T1D in prior work in DAISY (72). However, dairy intake contributes modestly to OCFA levels. These lipids primarily originate endogenously from adipocytes as well as from dietary intake of numerous foods including dairy, poultry, and fiber (70, 73, 74). Additional work is needed to understand connections between increasing methylation and increasing OCFA as well as the source of these lipid species.

A major strength of our study was the inclusion of DNAm measurements prior to T1D as well as the multi-omics work used to identify correlations between DNAm and gene expression as well metabolite levels. We measured DNAm before and after SV (ie, the appearance of IA) which builds on prior studies that have included DNAm measures after T1D and/or after IA onset only (14–16). A novel feature of our longitudinal methodology was our group*visit interaction modelling strategy that allowed us to identify changes in DNAm before and after the onset of IA, a critical window in T1D pathogenesis. These within individual effects are essential to understanding the etiology of T1D as they are robust to individual level confounders such as sex, genetic predisposition, and/or family history. Johnson et al (17) also used a longitudinal case-control analysis of T1D cases vs. unaffected controls in DAISY. In contrast, the current study design focused on individuals who developed IA and furthermore, tested for differences in DNAm post- vs pre-SV (group*visit interaction) rather than testing for differences in methylation by age (group*age interaction). Comparing the DMRs identified by this study versus Johnson et al (17), only two regions were located within 1 MB of each other–one on chr 6 ΔDMR 9 (28945322–28945493) in the current study vs chr 6 28973328-28973521 in Johnson et al (17), and one on chr 20 ΔDMR 2 (36148604–36149751) in the current study vs chr 20 36148954-36149232 in Johnson et al (17). Consistent with prior work, ΔDMR 9 and ΔDMR 2 were both associated with differential changes in DNAm in progressors relative to maintainers and/or reverters.

4.1 Limitations

We obtained DNAm from whole blood, which means we were unable to identify cell subtype specific effects. Similarly, our study focused on blood tissue only. DNAm changes within the blood may not reflect DNAm changes within other tissues that contribute to T1D, such as the pancreas. Due to advancements in technology during the study, DNAm was measured on two platforms. Individuals were randomly assigned to the platforms to minimize bias. We looked for cis-eQTMs. Given that it is possible that regions act over larger areas of the genome, we may have missed larger effects that occurred outside of our 500 KB window. Due to the small sample size, the eQTM was underpowered to identify FDR significant DMR vs gene transcript pairs. This limitation may explain lack of concordance between eQTM results and BIOS QTL results (μDMR4). Furthermore, among the two gene transcripts that were correlated with changes in methylation within ΔDMR1, gene expression data were only available at the post-SV visit. Therefore, it was not possible to determine whether gene expression also changed pre- versus post-SV. Finally, metabolite levels are influenced by age and dietary patterns. Although we adjusted for age, the large differences in age between the progressor group and the reverter and maintainer groups creates challenges in interpreting the metabolite vs methylation correlations. Additional work is need to replicate the metabolite vs DNAm regional effects.

5 Conclusion

T1D is an autoimmune disease characterized by immune mediated destruction of beta cells. Beta cell stress has been proposed as a mechanism connecting environmental perturbations such as infection, inflammation, diet, and increased insulin secretion to disease progression (75). Our EWAS identified DNAm candidates known to be modified by diabetes relevant environmental factors including diet and glucose levels (CLIP2, GNAS/GNAS-AS, MEG3). Our results also implicated genes (DLL1 and GNAS) with functional roles in both beta and immune cells. Our results build upon prior work by identifying specific areas of the genome where DNAm changes pre- and post-SV visits differentiated between reversion versus progression of IA. The correlation between changes in DNAm and changes in lipid levels reveal common connections between DMRs in different areas of the genome that may be related to disruptions in lipid metabolic pathways. Additional work is needed to replicate these findings, test for cell-specific changes in DNAm pre- vs post-seroconversion, and to identify modifiable factors that lead to these DNAm changes; ideally, the first step in the development of preventative strategies that delay or prevent progression of IA.

Data availability statement

The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found below: https://www.ncbi.nlm.nih.gov/, PRJNA597238.

Ethics statement

The studies involving humans were approved by Colorado Multi-institutional Review Board. The studies were conducted in accordance with the local legislation and institutional requirements. Written informed consent for participation in this study was provided by the participants’ legal guardians/next of kin.

Author contributions

PC: Writing – review & editing, Writing – original draft, Project administration, Methodology, Formal analysis, Conceptualization. LV: Writing – review & editing, Software, Methodology, Data curation. RJ: Writing – review & editing, Methodology, Data curation. TB: Writing – review & editing. AS: Writing – review & editing, Supervision. KK: Writing – review & editing, Supervision, Methodology. IY: Writing – review & editing, Supervision. TF: Writing – review & editing, Supervision. OF: Writing – review & editing, Resources, Data curation. MR: Writing – review & editing, Supervision, Resources, Funding acquisition. JN: Writing – review & editing, Supervision, Methodology, Funding acquisition, Conceptualization.

Funding

The author(s) declare financial support was received for the research, authorship, and/or publication of this article. This work was supported by National Institutes R01-DK104351, R01-DK32493, R21-AI142483, and P30-DK116073.

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

The author(s) declared that they were an editorial board member of Frontiers, at the time of submission. This had no impact on the peer review process and the final decision.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fimmu.2024.1345494/full#supplementary-material

References

1. Ziegler AG, Rewers M, Simell O, Simell T, Lempainen J, Steck A, et al. Seroconversion to multiple islet autoantibodies and risk of progression to diabetes in children. JAMA. (2013) 309:2473–9. doi: 10.1001/jama.2013.6285

PubMed Abstract | CrossRef Full Text | Google Scholar

2. Vehik K, Lynch KF, Schatz DA, Akolkar B, Hagopian W, Rewers M, et al. Reversion of beta-cell autoimmunity changes risk of type 1 diabetes: TEDDY study. Diabetes Care. (2016) 39:1535–42. doi: 10.2337/dc16-0181

PubMed Abstract | CrossRef Full Text | Google Scholar

3. Spencer KM, Tarn A, Dean BM, Lister J, Bottazzo GF. Fluctuating islet-cell autoimmunity in unaffected relatives of patients with insulin-dependent diabetes. Lancet. (1984) 1:764–6. doi: 10.1016/S0140-6736(84)91278-9

PubMed Abstract | CrossRef Full Text | Google Scholar

4. Yu J, Yu L, Bugawan TL, Erlich HA, Barriga K, Hoffman M, et al. Transient antiislet autoantibodies: infrequent occurrence and lack of association with "genetic" risk factors. J Clin Endocrinol Metab. (2000) 85:2421–8. doi: 10.1210/jc.85.7.2421

PubMed Abstract | CrossRef Full Text | Google Scholar

5. Kulmala P, Rahko J, Savola K, Vahasalo P, Veijola R, Sjoroos M, et al. Stability of autoantibodies and their relation to genetic and metabolic markers of Type I diabetes in initially unaffected schoolchildren. Diabetologia. (2000) 43:457–64. doi: 10.1007/s001250051329

PubMed Abstract | CrossRef Full Text | Google Scholar

6. Kimpimaki T, Kulmala P, Savola K, Kupila A, Korhonen S, Simell T, et al. Natural history of beta-cell autoimmunity in young children with increased genetic susceptibility to type 1 diabetes recruited from the general population. J Clin Endocrinol Metab. (2002) 87:4572–9. doi: 10.1210/jc.2002-020018

PubMed Abstract | CrossRef Full Text | Google Scholar

7. Robertson CC, Inshaw JRJ, Onengut-Gumuscu S, Chen WM, Santa Cruz DF, Yang H, et al. Fine-mapping, trans-ancestral and genomic analyses identify causal variants, cells, genes and drug targets for type 1 diabetes. Nat Genet. (2021) 53:962–71. doi: 10.1038/s41588-021-00880-5

PubMed Abstract | CrossRef Full Text | Google Scholar

8. Jerram ST, Leslie RD. The genetic architecture of type 1 diabetes. Genes (Basel). (2017) 8. doi: 10.3390/genes8080209

PubMed Abstract | CrossRef Full Text | Google Scholar

9. Patterson CC, Dahlquist GG, Gyurus E, Green A, Soltesz G, Group ES. Incidence trends for childhood type 1 diabetes in Europe during 1989-2003 and predicted new cases 2005-20: a multicentre prospective registration study. Lancet. (2009) 373:2027–33. doi: 10.1016/S0140-6736(09)60568-7

PubMed Abstract | CrossRef Full Text | Google Scholar

10. Ziegler AG, Pflueger M, Winkler C, Achenbach P, Akolkar B, Krischer JP, et al. Accelerated progression from islet autoimmunity to diabetes is causing the escalating incidence of type 1 diabetes in young children. J Autoimmun. (2011) 37:3–7. doi: 10.1016/j.jaut.2011.02.004

PubMed Abstract | CrossRef Full Text | Google Scholar

11. Lappalainen T, Greally JM. Associating cellular epigenetic models with human phenotypes. Nat Rev Genet. (2017) 18:441–51. doi: 10.1038/nrg.2017.32

PubMed Abstract | CrossRef Full Text | Google Scholar

12. Sawalha AH. Epigenetics and T-cell immunity. Autoimmunity. (2008) 41:245–52. doi: 10.1080/08916930802024145

PubMed Abstract | CrossRef Full Text | Google Scholar

13. Allan RS, Nutt SL. Deciphering the epigenetic code of T lymphocytes. Immunol Rev. (2014) 261:50–61. doi: 10.1111/imr.12207

PubMed Abstract | CrossRef Full Text | Google Scholar

14. Rakyan VK, Beyan H, Down TA, Hawa MI, Maslau S, Aden D, et al. Identification of type 1 diabetes-associated DNA methylation variable positions that precede disease diagnosis. PloS Genet. (2011) 7:e1002300. doi: 10.1371/journal.pgen.1002300

PubMed Abstract | CrossRef Full Text | Google Scholar

15. Paul DS, Teschendorff AE, Dang MA, Lowe R, Hawa MI, Ecker S, et al. Increased DNA methylation variability in type 1 diabetes across three immune effector cell types. Nat Commun. (2016) 7:13555. doi: 10.1038/ncomms13555

PubMed Abstract | CrossRef Full Text | Google Scholar

16. Elboudwarej E, Cole M, Briggs FB, Fouts A, Fain PR, Quach H, et al. Hypomethylation within gene promoter regions and type 1 diabetes in discordant monozygotic twins. J Autoimmun. (2016) 68:23–9. doi: 10.1016/j.jaut.2015.12.003

PubMed Abstract | CrossRef Full Text | Google Scholar

17. Johnson RK, Vanderlinden LA, Dong F, Carry PM, Seifert J, Waugh K, et al. Longitudinal DNA methylation differences precede type 1 diabetes. Sci Rep. (2020) 10:3721. doi: 10.1038/s41598-020-60758-0

PubMed Abstract | CrossRef Full Text | Google Scholar

18. Yu L, Robles DT, Abiru N, Kaur P, Rewers M, Kelemen K, et al. Early expression of antiinsulin autoantibodies of humans and the NOD mouse: evidence for early determination of subsequent diabetes. Proc Natl Acad Sci USA. (2000) 97:1701–6. doi: 10.1073/pnas.040556697

PubMed Abstract | CrossRef Full Text | Google Scholar

19. Yu L, Rewers M, Gianani R, Kawasaki E, Zhang Y, Verge C, et al. Antiislet autoantibodies usually develop sequentially rather than simultaneously. J Clin Endocrinol Metab. (1996) 81:4264–7. doi: 10.1210/jcem.81.12.8954025

PubMed Abstract | CrossRef Full Text | Google Scholar

20. Bonifacio E, Yu L, Williams AK, Eisenbarth GS, Bingley PJ, Marcovina SM, et al. Harmonization of glutamic acid decarboxylase and islet antigen-2 autoantibody assays for national institute of diabetes and digestive and kidney diseases consortia. J Clin Endocrinol Metab. (2010) 95:3360–7. doi: 10.1210/jc.2010-0293

PubMed Abstract | CrossRef Full Text | Google Scholar

21. Wenzlau JM, Juhl K, Yu L, Moua O, Sarkar SA, Gottlieb P, et al. The cation efflux transporter ZnT8 (Slc30A8) is a major autoantigen in human type 1 diabetes. Proc Natl Acad Sci USA. (2007) 104:17040–5. doi: 10.1073/pnas.0705894104

PubMed Abstract | CrossRef Full Text | Google Scholar

22. Aryee MJ, Jaffe AE, Corrada-Bravo H, Ladd-Acosta C, Feinberg AP, Hansen KD, et al. Minfi: a flexible and comprehensive Bioconductor package for the analysis of Infinium DNA methylation microarrays. Bioinformatics. (2014) 30:1363–9. doi: 10.1093/bioinformatics/btu049

PubMed Abstract | CrossRef Full Text | Google Scholar

23. Vanderlinden LA, Johnson RK, Carry PM, Dong F, DeMeo DL, Yang IV, et al. An effective processing pipeline for harmonizing DNA methylation data from Illumina's 450K and EPIC platforms for epidemiological studies. BMC Res Notes. (2021) 14:352. doi: 10.1186/s13104-021-05741-2

PubMed Abstract | CrossRef Full Text | Google Scholar

24. Carry PM, Waugh K, Vanderlinden LA, Johnson RK, Buckner T, Rewers M, et al. Changes in the co-expression of innate immunity genes during persistent islet autoimmunity are associated with progression of islet autoimmunity: the diabetes autoimmunity study in the young (DAISY). Diabetes. (2022) 71(9):2048–57. doi: 10.2337/figshare.20060480

PubMed Abstract | CrossRef Full Text | Google Scholar

25. Li B, Dewey CN. RSEM: accurate transcript quantification from RNA-Seq data with or without a reference genome. BMC Bioinf. (2011) 12:323. doi: 10.1186/1471-2105-12-323

CrossRef Full Text | Google Scholar

26. Love MI, Huber W, Anders S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol. (2014) 15:550. doi: 10.1186/s13059-014-0550-8

PubMed Abstract | CrossRef Full Text | Google Scholar

27. Risso D, Ngai J, Speed TP, Dudoit S. Normalization of RNA-seq data using factor analysis of control genes or samples. Nat Biotechnol. (2014) 32:896–902. doi: 10.1038/nbt.2931

PubMed Abstract | CrossRef Full Text | Google Scholar

28. Carry PM, Vanderlinden LA, Johnson RK, Buckner T, Fiehn O, Steck AK, et al. Phospholipid levels at seroconversion are associated with resolution of persistent islet autoimmunity: the diabetes autoimmunity study in the young. Diabetes. (2021) 70:1592–601. doi: 10.2337/db20-1251

PubMed Abstract | CrossRef Full Text | Google Scholar

29. Showalter MR, Nonnecke EB, Linderholm AL, et al. Obesogenic diets alter metabolism in mice. PLoS One. (2018) 13:e0190632.

PubMed Abstract | Google Scholar

30. Fiehn O, Wohlgemuth G, Scholz M, et al. Quality control for plant metabolomics: reporting MSI-compliant studies. Plant J. (2008) 53:691–704.

PubMed Abstract | Google Scholar

31. Cajka T, Smilowitz JT, Fiehn O. Validating quantitative untargeted lipidomics across nine liquid chromatography-high-resolution mass spectrometry platforms. Anal Chem. (2017) 89:12360–80.

PubMed Abstract | Google Scholar

32. Skogerson K, Wohlgemuth G, Barupal DK, Fiehn O. The volatile compound BinBase mass spectral database. BMC Bioinformatics. (2011) 12:321.

PubMed Abstract | Google Scholar

33. Tsugawa H, Cajka T, Kind T, et al. MS-DIAL: data-independent MS/MS deconvolution for comprehensive metabolome analysis. Nat Methods. (2015) 12:523–6.

PubMed Abstract | Google Scholar

34. Kind T, Liu KH, Lee DY, DeFelice B, Meissen JK, Fiehn O. LipidBlast in silico tandem mass spectrometry database for lipid identification. Nat Methods. (2013) 10:755–8.

PubMed Abstract | Google Scholar

35. Fan S, Kind T, Cajka T, et al. Systematic error removal using random forest for normalizing large-scale untargeted lipidomics data. Anal Chem. (2019) 91:3590–6.

PubMed Abstract | Google Scholar

36. Buckner T, Johnson RK, Vanderlinden LA, Carry PM, Romero A, Onengut-Gumuscu S, et al. Genome-wide analysis of oxylipins and oxylipin profiles in a pediatric population. Front Nutr. (2023) 10:1040993. doi: 10.3389/fnut.2023.1040993

PubMed Abstract | CrossRef Full Text | Google Scholar

37. Benjamini Y, Hochberg Y. Controlling the false discovery rate: a practical and powerful approach to multiple testing. J R Stat Soc Ser B (Methodological). (1995) 57:289–300. doi: 10.1111/j.2517-6161.1995.tb02031.x

CrossRef Full Text | Google Scholar

38. Pedersen BS, Schwartz DA, Yang IV, Kechris KJ. Comb-p: software for combining, analyzing, grouping and correcting spatially correlated P-values. Bioinformatics. (2012) 28:2986–8. doi: 10.1093/bioinformatics/bts545

PubMed Abstract | CrossRef Full Text | Google Scholar

39. Dreja T, Jovanovic Z, Rasche A, Kluge R, Herwig R, Tung YC, et al. Diet-induced gene expression of isolated pancreatic islets from a polygenic mouse model of the metabolic syndrome. Diabetologia. (2010) 53:309–20. doi: 10.1007/s00125-009-1576-4

PubMed Abstract | CrossRef Full Text | Google Scholar

40. Diabetes Genetics Initiative of Broad Institute of H, Mit LU, Novartis Institutes of BioMedical R, Saxena R, Voight BF, Lyssenko V, et al. Genome-wide association analysis identifies loci for type 2 diabetes and triglyceride levels. Science. (2007) 316:1331–6. doi: 10.1126/science.1142358

PubMed Abstract | CrossRef Full Text | Google Scholar

41. Zeggini E, Weedon MN, Lindgren CM, Frayling TM, Elliott KS, Lango H, et al. Replication of genome-wide association signals in UK samples reveals risk loci for type 2 diabetes. Science. (2007) 316:1336–41. doi: 10.1126/science.1142364

PubMed Abstract | CrossRef Full Text | Google Scholar

42. Consortium EP. An integrated encyclopedia of DNA elements in the human genome. Nature. (2012) 489:57–74. doi: 10.1038/nature11247

PubMed Abstract | CrossRef Full Text | Google Scholar

43. Alliance of Genome Resources C. Alliance of Genome Resources Portal: unified model organism research platform. Nucleic Acids Res. (2020) 48:D650–D8. doi: 10.1093/nar/gkz813

PubMed Abstract | CrossRef Full Text | Google Scholar

44. Kim W, Shin YK, Kim BJ, Egan JM. Notch signaling in pancreatic endocrine cell and diabetes. Biochem Biophys Res Commun. (2010) 392:247–51. doi: 10.1016/j.bbrc.2009.12.115

PubMed Abstract | CrossRef Full Text | Google Scholar

45. Garis M, Garrett-Sinha LA. Notch signaling in B cell immune responses. Front Immunol. (2020) 11:609324. doi: 10.3389/fimmu.2020.609324

PubMed Abstract | CrossRef Full Text | Google Scholar

46. Apelqvist A, Li H, Sommer L, Beatus P, Anderson DJ, Honjo T, et al. Notch signalling controls pancreatic cell differentiation. Nature. (1999) 400:877–81. doi: 10.1038/23716

PubMed Abstract | CrossRef Full Text | Google Scholar

47. Rubey M, Chhabra NF, Gradinger D, Sanz-Moreno A, Lickert H, Przemeck GKH, et al. DLL1- and DLL4-mediated notch signaling is essential for adult pancreatic islet homeostasis. Diabetes. (2020) 69:915–26. doi: 10.2337/db19-0795

PubMed Abstract | CrossRef Full Text | Google Scholar

48. Dror V, Nguyen V, Walia P, Kalynyak TB, Hill JA, Johnson JD. Notch signalling suppresses apoptosis in adult human and mouse pancreatic islet cells. Diabetologia. (2007) 50:2504–15. doi: 10.1007/s00125-007-0835-5

PubMed Abstract | CrossRef Full Text | Google Scholar

49. Al-Shami A, Crisostomo J, Wilkins C, Xu N, Humphries J, Chang WC, et al. Integrin-alpha FG-GAP repeat-containing protein 2 is critical for normal B cell differentiation and controls disease development in a lupus model. J Immunol. (2013) 191:3789–98. doi: 10.4049/jimmunol.1203534

PubMed Abstract | CrossRef Full Text | Google Scholar

50. Bonder MJ, Luijk R, Zhernakova DV, Moed M, Deelen P, Vermaat M, et al. Disease variants alter transcription factor levels and methylation of their binding sites. Nat Genet. (2017) 49:131–8. doi: 10.1038/ng.3721

PubMed Abstract | CrossRef Full Text | Google Scholar

51. Chen D, Zhang A, Fang M, Fang R, Ge J, Jiang Y, et al. Increased methylation at differentially methylated region of GNAS in infants born to gestational diabetes. BMC Med Genet. (2014) 15:108. doi: 10.1186/s12881-014-0108-3

PubMed Abstract | CrossRef Full Text | Google Scholar

52. Lumey LH, Stein AD, Kahn HS, van der Pal-de Bruin KM, Blauw GJ, Zybert PA, et al. Cohort profile: the Dutch Hunger Winter families study. Int J Epidemiol. (2007) 36:1196–204. doi: 10.1093/ije/dym126

PubMed Abstract | CrossRef Full Text | Google Scholar

53. Tobi EW, Lumey LH, Talens RP, Kremer D, Putter H, Stein AD, et al. DNA methylation differences after exposure to prenatal famine are common and timing- and sex-specific. Hum Mol Genet. (2009) 18:4046–53. doi: 10.1093/hmg/ddp353

PubMed Abstract | CrossRef Full Text | Google Scholar

54. Weinstein LS, Xie T, Zhang QH, Chen M. Studies of the regulation and function of the Gs alpha gene Gnas using gene targeting technology. Pharmacol Ther. (2007) 115:271–91. doi: 10.1016/j.pharmthera.2007.03.013

PubMed Abstract | CrossRef Full Text | Google Scholar

55. Miyoshi N, Wagatsuma H, Wakana S, Shiroishi T, Nomura M, Aisaka K, et al. Identification of an imprinted gene, Meg3/Gtl2 and its human homologue MEG3, first mapped on mouse distal chromosome 12 and human chromosome 14q. Genes Cells. (2000) 5:211–20. doi: 10.1046/j.1365-2443.2000.00320.x

PubMed Abstract | CrossRef Full Text | Google Scholar

56. Taneera J, Dhaiban S, Mohammed AK, Mukhopadhyay D, Aljaibeji H, Sulaiman N, et al. GNAS gene is an important regulator of insulin secretory capacity in pancreatic beta-cells. Gene. (2019) 715:144028. doi: 10.1016/j.gene.2019.144028

PubMed Abstract | CrossRef Full Text | Google Scholar

57. Lee J, Kim TH, Murray F, Li X, Choi SS, Broide DH, et al. Cyclic AMP concentrations in dendritic cells induce and regulate Th2 immunity and allergic asthma. Proc Natl Acad Sci USA. (2015) 112:1529–34. doi: 10.1073/pnas.1417972112

PubMed Abstract | CrossRef Full Text | Google Scholar

58. Shah A. The pathologic and clinical intersection of atopic and autoimmune disease. Curr Allergy Asthma Rep. (2012) 12:520–9. doi: 10.1007/s11882-012-0293-0

PubMed Abstract | CrossRef Full Text | Google Scholar

59. Smew AI, Lundholm C, Savendahl L, Lichtenstein P, Almqvist C. Familial coaggregation of asthma and type 1 diabetes in children. JAMA Netw Open. (2020) 3:e200834. doi: 10.1001/jamanetworkopen.2020.0834

PubMed Abstract | CrossRef Full Text | Google Scholar

60. Metsala J, Lundqvist A, Virta LJ, Kaila M, Gissler M, Virtanen SM, et al. The association between asthma and type 1 diabetes: a paediatric case-cohort study in Finland, years 1981-2009. Int J Epidemiol. (2018) 47:409–16. doi: 10.1093/ije/dyx245

PubMed Abstract | CrossRef Full Text | Google Scholar

61. Yun HD, Knoebel E, Fenta Y, Gabriel SE, Leibson CL, Loftus EV Jr., et al. Asthma and proinflammatory conditions: a population-based retrospective matched cohort study. Mayo Clin Proc. (2012) 87:953–60. doi: 10.1016/j.mayocp.2012.05.020

PubMed Abstract | CrossRef Full Text | Google Scholar

62. La Torre D, Seppanen-Laakso T, Larsson HE, Hyotylainen T, Ivarsson SA, Lernmark A, et al. Decreased cord-blood phospholipids in young age-at-onset type 1 diabetes. Diabetes. (2013) 62:3951–6. doi: 10.2337/db13-0215

PubMed Abstract | CrossRef Full Text | Google Scholar

63. Lamichhane S, Ahonen L, Dyrlund TS, Kemppainen E, Siljander H, Hyoty H, et al. Dynamics of plasma lipidome in progression to islet autoimmunity and type 1 diabetes - type 1 diabetes prediction and prevention study (DIPP). Sci Rep. (2018) 8:10635. doi: 10.1038/s41598-018-28907-8

PubMed Abstract | CrossRef Full Text | Google Scholar

64. Oresic M, Gopalacharyulu P, Mykkanen J, Lietzen N, Makinen M, Nygren H, et al. Cord serum lipidome in prediction of islet autoimmunity and type 1 diabetes. Diabetes. (2013) 62:3268–74. doi: 10.2337/db13-0159

PubMed Abstract | CrossRef Full Text | Google Scholar

65. Li Q, Parikh H, Butterworth MD, Lernmark A, Hagopian W, Rewers M, et al. Longitudinal metabolome-wide signals prior to the appearance of a first islet autoantibody in children participating in the TEDDY study. Diabetes. (2020) 69:465–76. doi: 10.2337/db19-0756

PubMed Abstract | CrossRef Full Text | Google Scholar

66. Johnson RK, Vanderlinden L, DeFelice BC, Kechris K, Uusitalo U, Fiehn O, et al. Metabolite-related dietary patterns and the development of islet autoimmunity. Sci Rep. (2019) 9:14819.

PubMed Abstract | Google Scholar

67. Oresic M, Simell S, Sysi-Aho M, Nanto-Salonen K, Seppanen-Laakso T, Parikka V, et al. Dysregulation of lipid and amino acid metabolism precedes islet autoimmunity in children who later progress to type 1 diabetes. J Exp Med. (2008) 205:2975–84. doi: 10.1084/jem.20081800

PubMed Abstract | CrossRef Full Text | Google Scholar

68. Lamichhane S, Ahonen L, Dyrlund TS, Dickens AM, Siljander H, Hyoty H, et al. Cord-blood lipidome in progression to islet autoimmunity and type 1 diabetes. Biomolecules. (2019) 9. doi: 10.3390/biom9010033

PubMed Abstract | CrossRef Full Text | Google Scholar

69. Huynh K, Barlow CK, Jayawardana KS, Weir JM, Mellett NA, Cinel M, et al. High-throughput plasma lipidomics: detailed mapping of the associations with cardiometabolic risk factors. Cell Chem Biol. (2019) 26:71–84 e4. doi: 10.1016/j.chembiol.2018.10.008

PubMed Abstract | CrossRef Full Text | Google Scholar

70. Prada M, Wittenbecher C, Eichelmann F, Wernitz A, Drouin-Chartier JP, Schulze MB. Association of the odd-chain fatty acid content in lipid groups with type 2 diabetes risk: A targeted analysis of lipidomics data in the EPIC-Potsdam cohort. Clin Nutr. (2021) 40:4988–99. doi: 10.1016/j.clnu.2021.06.006

PubMed Abstract | CrossRef Full Text | Google Scholar

71. Pflueger M, Seppanen-Laakso T, Suortti T, Hyotylainen T, Achenbach P, Bonifacio E, et al. Age- and islet autoimmunity-associated differences in amino acid and lipid metabolites in children at risk for type 1 diabetes. Diabetes. (2011) 60:2740–7. doi: 10.2337/db10-1652

PubMed Abstract | CrossRef Full Text | Google Scholar

72. Lamb MM, Miller M, Seifert JA, Frederiksen B, Kroehl M, Rewers M, et al. The effect of childhood cow's milk intake and HLA-DR genotype on risk of islet autoimmunity and type 1 diabetes: the Diabetes Autoimmunity Study in the Young. Pediatr Diabetes. (2015) 16:31–8. doi: 10.1111/pedi.2015.16.issue-1

PubMed Abstract | CrossRef Full Text | Google Scholar

73. Green CR, Wallace M, Divakaruni AS, Phillips SA, Murphy AN, Ciaraldi TP, et al. Branched-chain amino acid catabolism fuels adipocyte differentiation and lipogenesis. Nat Chem Biol. (2016) 12:15–21. doi: 10.1038/nchembio.1961

PubMed Abstract | CrossRef Full Text | Google Scholar

74. Dąbrowski G, Konopka I. Update on food sources and biological activity of odd-chain, branched and cyclic fatty acids—-A review. Trends Food Sci Technol. (2021) 119:514–29. doi: 10.1016/j.tifs.2021.12.019

CrossRef Full Text | Google Scholar

75. Rewers M, Ludvigsson J. Environmental risk factors for type 1 diabetes. Lancet. (2016) 387:2340–8. doi: 10.1016/S0140-6736(16)30507-4

PubMed Abstract | CrossRef Full Text | Google Scholar