Dynamics of Plasma Lipidome in Progression to Islet Autoimmunity and Type 1 Diabetes – Type 1 Diabetes Prediction and Prevention Study (DIPP)

Type 1 diabetes (T1D) is one of the most prevalent autoimmune diseases among children in Western countries. Earlier metabolomics studies suggest that T1D is preceded by dysregulation of lipid metabolism. Here we used a lipidomics approach to analyze molecular lipids in a prospective series of 428 plasma samples from 40 children who progressed to T1D (PT1D), 40 children who developed at least a single islet autoantibody but did not progress to T1D during the follow-up (P1Ab) and 40 matched controls (CTR). Sphingomyelins were found to be persistently downregulated in PT1D when compared to the P1Ab and CTR groups. Triacylglycerols and phosphatidylcholines were mainly downregulated in PT1D as compared to P1Ab at the age of 3 months. Our study suggests that distinct lipidomic signatures characterize children who progressed to islet autoimmunity or overt T1D, which may be helpful in the identification of at-risk children before the initiation of autoimmunity.


Results
Age related changes in the lipidome. We analyzed plasma lipidomics prospectively in three study groups: 40 progressors to T1D (PT1D), 40 who progressed to at least a single Ab but not to T1D during the follow-up (P1Ab), and 40 controls (CTR) subjects who remained islet autoantibody negative during the follow-up until the age of 15. For each child, up to six time points were included in the analysis, corresponding to the ages of 3, 6, 12, 18, 24 or 36 months (Fig. 1). None of the children were diagnosed with T1D at the time of sample withdrawal. The lipidomics dataset (n = 428) used for the analysis included identified lipids from the following lipid classes: cholesterol esters (CE), diacylglycerols (DG), lysophosphatidylcholines (LPC), phosphatidylcholines (PC), phosphatidylethanolamines (PE), sphingomyelins (SM) and triacylglycerols (TG).
Principal Components Analysis (PCA) 13 of the preprocessed lipidomics data revealed a clear age-related trend. In order to specifically determine the age-related associations in the plasma lipidome, analysis of variance (ANOVA)-simultaneous component analysis (ASCA) 14 was performed with factors age, gender and case status (PT1D, P1Ab, CTR) including their interactions separating the variation into contributions from the different factors. Age had the strongest effect (6.01, p = 0.0010) on the plasma lipidomic profile when compared to the case group (0.89, p = 0.0050) and gender (0.60, p = 0.0020). However, the interactions between these factors remained insignificant. The score plot in Fig. 2a highlights the age-related trend and the corresponding loadings plot (Fig. 2b) explains the class-specific lipid alteration related to age. Specifically, the loadings depicts that the levels of SMs and LPCs were elevated at the early time points (3, 6, 12 months) while CEs and PEs appeared to be higher at later time points (18, 24, 36 months). In contrast, the levels of TGs and PCs varied over time with no clear age-related patterns. Figure 1. An overview of the study setting. The plasma samples for lipidomic analysis were obtained from the Finnish Type 1 Diabetes Prevention and Prediction Study (DIPP). The study cohort comprises the samples from three groups: children who progressed to T1D during the follow-up (PT1D), who developed at least a single islet autoantibody but did not progress to T1D during the follow-up (P1Ab), and control (CTR) subjects who remained islet autoantibody negative during the follow-up until the age of 15 years. The three study groups were matched by HLA-associated diabetes risk, gender and period of birth. For each child, longitudinal samples were obtained corresponding to the ages of 3, 6, 12, 18, 24, and 36 months. These age groups were selected with the objective of understanding the dynamics of the lipidome preceding the overt T1D. For each age cohort and study group, number of autoantibody positive children is marked.
SCIenTIfIC REPORTS | (2018) 8:10635 | DOI: 10.1038/s41598-018-28907-8 Plasma lipidome in development of autoimmunity and progression to T1D. In order to minimize the confounding effect of age, data were analyzed separately in different age cohorts. Partial least squares discriminant analysis (PLS-DA) 15 was performed to examine the differences between the three sample groups using combinations of various categorical variables (e.g., PT1D vs. P1Ab, PT1D vs. CTR) in each age cohort. Figure 3 shows the VIP (variable importance for projection) vs. regression coefficient plot as obtained from the PLS-DA analysis, between PT1D (n = 15) and CTR (n = 18) at 3 months of age (AUROC = 0.78). The result show lower levels of SMs and higher levels of some PCs in children who later progressed to T1D as compared to CTR. Among the PCs, particularly those with a higher number of double bonds (>1) showed an increase in PT1D as compared to CTR ( Supplementary Fig. 1). In addition, CE (20:4) and PE (36:4) appeared to be higher in PT1D, but no clear patterns were seen for the TGs and LPCs. For example, LPC (16:1) was elevated while LPC (16:0e) was downregulated in PT1D. Based on the VIP scores, SM (d36:2), SM (d36:1), SM (d34:1), SM (d38:2) and PC (40:4), PC (O-34:2), PC (34:2) and PC (O-36:2) were top ranked discriminating lipids in the supervised analysis. Univariate analysis showed that at the age of 3 months, a total of 12 lipids were altered between PT1D and CTR according to the nominal p-value < 0.05. These lipids included four PCs, seven SMs and one TG (Supplementary  Table 1). However, none of these lipids passed significance at the selected false discovery rate (FDR) threshold of 0.1. The adjusted p-values for SM (d36:1), SM (d34:1), and SM (d36:2) were 0.15, 0.15, and 0.20, respectively. Here each sample is represented by a point and colored according to the age (red diamond: 3 month of age, green square: 6 month of age, cyan triangle down: 12 months of age, light yellow circle: 18 months, blue square: 24 months, grey white circle 36 months). These ages of the participants are also marked on the x-axis while y-axis represents the sample score. Samples with similar score are clustered together. (b) The corresponding PC1 loading plot. The loadings explain the pattern seen in the score plot which provides the means to interpret the class specific lipid alteration related to age. X-axis is the variable order (i.e. lipids) and y-axis represents the lipids pattern corresponding to the score plot. Each color in the loadings represents a lipid class: cholesterol ester (CE), diacylglycerol (DG), lysophosphatidylcholine (LPC), phosphatidylcholine (PC), phosphatidylethanolamine (PE), sphingomyelin (SM) and triacylglycerol (TG). The loading depicts that the levels of SMs (green circle's) and LPCs (light blue circle's) were elevated at early time points (3, 6, 12 months) while CEs and PEs appeared to be higher in later time point (18, 24, 36 months PLS-DA analysis comparing PT1D (n = 15) and P1Ab (n = 18) at the age of 3 months revealed higher levels of SMs and TGs in the P1Ab group (AUROC = 0.67; Supplementary Fig. 2a). Univariate analysis showed that a total of 45 lipids being different between PT1D and P1Ab (nominal p-value < 0.05). These lipids included one LPC, 12 PCs, one PE, nine SMs and 22 TGs. All of these lipids passed the FDR threshold of 0.1 (  (Fig. 4a). However, no clear pattern with respect to the acyl chain composition was observed in any specific lipid class.
Next, PLS-DA analysis between P1Ab (n = 18) and CTR (n = 18) at 3 months of age revealed increased levels of PCs and PEs in the P1Ab group (AUROC = 0.65; Supplementary Fig. 2b). However, no clear patterns were observed for the other lipid classes. Analysis of individual lipids revealed that a total of 45 lipids were altered between the two study groups (nominal p-value < 0.05). These lipids included two CEs, 38 PCs, four PEs, and two TGs. Fourty-four of these lipids passed the selected FDR threshold of 0.1 (Supplementary Table 3). The heatmap highlights that most of the PCs, with the exception of PC (33:0) were increased in the P1Ab group. ASCA analysis also revealed that the study group had a significant effect (p = 0.022) on the plasma lipidome profile of the children while the effect of gender (p = 0.60) and its interaction (p = 0.95) remained insignificant in the 3-month cohort.
Similarly, we analyzed the lipid concentration differences between the three study groups in the remaining longitudinal series at 6, 12, 18, 24 and 36 months. There was no persistent trend with respect of lipid differences between the groups (Fig. 4a), with the exception for SMs which were persistently low in PT1D as compared to P1Ab and CTR (Fig. 4b). Low levels of SMs were particularly observed at 6, 12 and 18 months (Fig. 4b). However, only one of the SMs, SM (d36:2), differed significantly (adjusted p < 0.1) between the PT1D and CTR at 6 and 18 months of age (Fig. 4c). In addition, PC (33:0), TG (56:5) and TG (18:2/18:1/18:1) passed the FDR threshold of 0.1 at 18 months of age (Fig. 4a). Furthermore, PLS-DA analysis performed between the three study groups (PT1D vs. CTR, PT1D vs. P1Ab, CTR vs. P1Ab) in the remaining longitudinal series at the age of 6, 12, 18, 24 and 36 months also showed poor classification accuracy (based on AUROC value, i.e. less than 50% classified accurately).
Lipidome in progression to type 1 diabetes. Next, we sought to determine the total lipid class specific changes in the longitudinal series of samples. Individual lipid concentrations in each lipid class were added, then subsequent data analysis considering each lipid class as a variable was performed. We found clear alterations in phospholipids, in particular SMs were significantly reduced in the PT1D group at the ages of 3 and 18 months when compared to the P1Ab group (Fig. 5a) and at the age of 3 month when compared to CTR (Fig. 5b). TGs and PCs were significantly lower in PT1D than in P1Ab at the age of 3 months (Fig. 5a,c). PEs and PCs were also higher in P1Ab than in CTR at the age of 3 months.
Lipidome changes before the emergence of islet autoantibodies. In order to examine whether the circulating lipids were associated with the appearance of islet autoantibodies, the plasma lipidome in PT1D with one Ab (n = 21) and P1Ab (n = 27) groups were compared between the samples obtained before and after the appearance of the first islet autoantibody, respectively. We applied the multi-level approach to exploit the paired data structure. Multi-level analysis revealed several differences before and after the appearance of islet autoantibodies. Specifically, SMs (SM (d36:1), SM (d36:2)) and LPCs (LPC (18:0), LPC (18:2), LPC (16:0)) were higher while CEs (CE (20:4), CE (20:5)) were lower before seroconversion than after it in both PT1D and P1Ab (Fig. 6).
As an exception to this pattern, LPC (20:3) appeared to be elevated after the emergence of the first autoantibody in P1Ab. Next, we set out to examine the individual lipid profiles by univariate analysis. The pairwise analysis revealed a total of 27 (two CEs, three LPCs, 11PCs, five SMs, six TGs) and 16 (two CEs, two LPCs, six PCs, four SMs, two TGs) lipids, respectively, being altered in both study groups (nominal p-value < 0.05, Supplementary  Tables 4 and 5). However, only CE (20:5) passed the level of significance at the selected FDR threshold of 0.1 in P1Ab.
We also studied whether lipidome profiles associated with the appearance of specific islet autoantibodies. The plasma lipidomes were compared before and after the emergence of first insulin autoantibodies (IAA) in P1Ab and PT1D with one Ab. Pairwise comparisons were also performed by combining other autoantibodies as a group (islet cell antibodies (ICA), islet antigen 2 autoantibodies (IA-2A), and GAD autoantibodies (GADA) in the remaining children. The pairwise comparison revealed a total of 17 (in PT1D with first IAA), one (in PT1D with all others autoantibodies excluding IAA), 22 (in P1Ab with first IAA), and 11 (in P1Ab with all others autoantibodies excluding IAA) lipids being altered at a nominal p-value < 0.05 before and after the appearance of specific  Tables 6-9). However, none of the lipid species passed the significance at the FDR threshold of 0.1.
Next, we analyzed the lipid concentration differences between the specific islet (IAA vs others including IA-2, GADA, ICA) autoantibodies study groups in the longitudinal series at 3, 6, 12, 18, 24 and 36 months. There were no significant (FDR threshold of 0.1) lipid differences between the specific islet groups, with the exception for 3 months of age, where we found a total of 23 lipids (two CEs, one PC, seven SMs, 13 TGs) being altered when comparing PT1D with IAA to P1Ab with IAA (Supplementary Table 10). All of these lipids passed significance at the FDR threshold of 0.1. SMs, TGs and PC were lower while CEs appeared to be higher at the age of 3 months in the PT1D group with IAA than in the P1Ab group with IAA as the first autoantibody ( Supplementary Fig. 3).

Discussion
Here we showed that specific lipid changes precede islet autoimmunity and T1D. Our results revealed that the children who progress to T1D in the follow-up tend to have a distinct and persistently dyregulated lipid profile as compared to those who later progress to islet autoimmunity but not to T1D. Our observations also show that plasma lipids, in particular CEs, PEs, SMs and LPCs, are strongly affected by age during the early development. This is in agreement with earlier studies which suggest that gender, age, and body mass index have significant impact on the longitudinal trajectories of plasma lipids 16,17 . Several studies depict that serum total cholesterol and TGs gradually increase with age during the early development 18,19 , but the data specifically from infants are limited to phospholipids such as SMs and LPCs 20 . Our findings are not unexpected however, because various host related and environmental factors including dietary changes during the development can impact the lipid profiles over time 21 .
We found that differences in the lipidomes between the three study groups were most pronounced at the earliest age of 3 months. SM levels were low inPT1D as compared to the other groups. At the level of individual lipid species, most significant and persistent differences were observed for SM (d36:2), which differed between PT1D and CTR, as well as between PT1D and P1Ab at 3 months of age. Most of the other SMs also tended to be persistently lower in PT1D than in both CTR and P1Ab. Interestingly, the lipid-related associations with progression to T1D appear to be the strongest in the group of children who first developed positivity for IAA, which are also the most common first autoantibody detected in young children who later progress to T1D 5 and tends to be associated with higher risk of progression as compared to other autoantibodies 22,23 .
Previous studies by us and others have also reported decreased serum SM levels in children who later progressed to T1D as well as in children with newly diagnosed T1D 24,25 . Sphingolipids are among the most abundant choline containing phospholipids in the mammalian cells, tissues and in the circulation. Sphingolipids have been shown to be potent regulators of immune cell activity 26 . In addition, sphingolipids are crucial for plasma The figure shows that SMs were downregulated in the PT1D group at the age of 3 months compared to P1Ab and CTR (adjusted p = 0.047). TGs were also downregulated in PT1D as compared to P1Ab (adjusted p = 0.079) at the age of 3 months. Furthermore, PEs and PCs were also higher in P1Ab than in CTR at the age of 3 months. There was no persistent trend with respect of total lipid concentration between CTR, P1Ab and PT1D in other the perspective series of samples.  26 . Our results suggest that SM metabolism is disrupted among the T1D progressors from early age, already preceding the appearance of islet autoantibodies.
An important consideration in the search of metabolic signature of T1D progression is that blood in young infants is sampled in non-fasting state. In the earlier study, SMs were found to be least variable within-individuals in nonfasted serum samples 12 . This suggests that SMs may also hold promise as predictive metabolic signature of T1D and thus complementing genetic testing and islet autoantibody determination. Clearly more and larger studies in heterogeneous populations are needed to examine the complementary diagnostic value of SMs.
The total TGs and PCs were downregulated in PT1D as compared to P1Ab at 3 months of age. Consistent with these observations, other studies have reported perturbations in phospholipids during the early infancy among the T1D progressors 10,24,25 as well as decreased TGs in newborn infants who later progressed to T1D 24 . Choline regulates secretion of triglyceride rich lipoprotein particles and its deficiency has been linked to lower serum TG concentration 28 . Since PCs are a major source of choline; an important epigenetic regulator in the body 29 , it is conceivable that alterations in the choline metabolism during infancy, e.g. by diet, may affect the PC levels, resulting in the decreased levels of TGs as well as PCs. In addition, the decreased levels of TGs may also suggest increased energy demand in PT1D group as compared to P1Ab 30 . TGs are concentrated stores of metabolic energy in biological tissues 31 . During energy demand, TGs are hydrolyzed into free fatty acids that are taken up by other tissues or cells to meet the energy requirements of the organism 32 .
The appearance of autoantibodies was associated with down-regulation of SMs and LPCs and up-regulation of CEs. LPC is the main phospholipid component of oxidized low density lipoprotein (oxLDL) produced as a hydrolysis product of PC 33 . LPC regulates various intracellular signaling events triggered by the oxLDL, including modulation of the toll like receptor (TLR) signaling pathway and activation of macrophages which are the hallmarks of the inflammatory and immunological processes 33,34 . These findings are in line with earlier observations that LPC dysregulation precedes the appearance of islet autoantibodies 24 .
Up-regulation of CEs following seroconversion to positivity for islet autoantibodies is a novel finding. Previous studies have provided evidence for the inflammatory effects of cholesterol accumulation 35 . Cholesterol accumulation activated the enzyme lecithin cholesterol acyltransferase (LCAT) that modifies cholesterol to CE. The esterified cholesterols are subsequently taken up by the liver and discharged through bile excretion 36 . An efficient cholesterol esterification would improve the reverse cholesterol transport leading to a reduction in inflammatory lesions 37 . In contrast, reduced cholesterol efflux may trigger pro-inflammatory signaling pathways such as TLR signaling, which results in the modulation of inflammatory responses 38 . However, understanding of the role of cholesterol homeostasis with respect to islet autoimmunity is currently limited.
Taken together, our findings support earlier studies by us and others, which suggest that dysregulation of lipid metabolism precedes islet autoimmunity and T1D, as well as add a new knowledge that the observed changes are specifically and persistently associated with progression to overt disease. The distinct lipidomics profile which is associated with the progression to T1D appears to be most pronounced in very young children, i.e. at 3 months of age in our study. This suggests (1) that diversification of dietary patterns already during the development may mask some of the T1D-associated metabolic signatures and (2) that metabolite-based prediction of T1D may be most feasible in early life, already prior to the appearance of islet autoantibodies. Already in young infants, the lipidomics signatures of children who progress to T1D later in life are clearly different from those that do seroconvert to a single islet autoantibody but do not progress to overt disease. Since several lipids found dysregulated in the PT1D group displayed an opposite pattern in the P1Ab group in relation to CTR, our findings suggest that PT1D and P1Ab have distinct metabotype. Since the current stratifications in studies aimed at T1D prevention are based on the detection of islet autoantibodies, lipidomics profiles may thus provide a valuable complementary tool for identifying children at highest risk of progression to T1D.

Methods
Sample and study design. The samples have been obtained from the Finnish Type 1 Diabetes Prevention and Prediction Study (DIPP) 39 . The DIPP study has screened more than 220,000 newborn infants for HLAconferred susceptibility to T1D in three university hospitals Turku, Tampere, and Oulu in Finland 40 . The subjects involved in the current study were chosen from the subset of DIPP children from the Tampere study center. Here, five longitudinal samples per child were collected between 1998 and 2012. For each child, longitudinal samples were obtained corresponding to either of the ages of 3, 6, 12, 18, 24, and 36 months. This cohort comprises the samples from 120 children: 40 progressors to T1D (PT1D), 40 who tested positive for at least one Ab in a minimum of two consecutive samples but did not progress to clinical T1D during the follow-up (P1Ab), and 40 controls (CTR) subjects who remained islet autoantibody negative during the follow-up until the age of 15. The three study groups were matched by HLA-associated diabetes risk, gender and period of birth. Selected characteristics of the study subjects are shown in Table 1. In total 428 non-fasting blood samples were collected for this study. Plasma was separated within 30 minutes after the blood collection by centrifugation at 1600 g for 20 minutes at room temperature. The plasma samples were stored at −80 °C until analyzed.
Detection of islet autoantibodies. The participants with HLA-conferred genetic susceptibility were prospectively observed for the appearance of T1D associated autoantibodies (islet cell antibodies (ICA), insulin autoantibodies (IAA), islet antigen 2 autoantibodies (IA-2A), and GAD autoantibodies (GADA). These    autoantibodies were measured in the Diabetes Research Laboratory, University of Oulu from the plasma samples taken at each follow-up visit as described 43 . ICA were detected with the use of indirect immunofluorescence, whereas the other three autoantibodies were quantified with the use of specific radio binding assays 44  The samples were analysed using an ultra-high-performance liquid chromatography quadrupole time-of-flight mass spectrometry method (UHPLC-Q-TOF-MS), which has been presented in detail previously 46 . Briefly, the UHPLC system used in this work was a 1290 Infinity system from Agilent Technologies (Santa Clara, CA, USA). The system was equipped with a multi sampler (maintained at 10 °C), a quaternary solvent manager and a column thermostat (maintained at 50 °C). Separations were performed on an ACQUITY UPLC ® BEH C18 column (2.1 mm × 100 mm, particle size 1.7 µm) by Waters (Milford, USA).
The mass spectrometer coupled to the UHPLC was a 6550 iFunnel quadrupole time of flight (Q-TOF) from Agilent Technologies interfaced with a dual jet stream electrospray (dual ESI) ion source. All analyses were performed in positive ion mode and MassHunter B.06.01 (Agilent Technologies) was used for all data acquisition. Quality control was performed throughout the dataset by including blanks, pure standard samples, extracted standard samples and control plasma samples. Relative standard deviations (%RSDs) for retention times and peak areas for lipid standards representing each lipid class in the control plasma samples (n = 37) and in the pure standard samples (n = 37) were calculated. The %RSDs for the retention times were on average below 0.2% for both the control plasma samples and for the pure standard samples. The peak areas were normalized to the peak area of PC (16:0/d30/18:1) in each showing %RSDs within accepted analytical limits at averages of 13% and 14% for the control plasma samples and for the pure standard samples, respectively. This shows that the method is reliable and repeatable throughout the sample set. Statistical analysis. All statistical analyses were based on log2-transformed intensity data. The transformed data were mean centered and auto scaled prior to multivariate analysis to improve the global interpretability. The multivariate analysis was done using the PLS Toolbox 8.2.1 (Eigenvector Research Inc., Manson, WA, USA) in MATLAB 2017b (Mathworks, Inc., Natick, MA, USA). PCA of the pre-processed data was initially performed to highlight patterns and to emphasize variation in the dataset. ANOVA-simultaneous component analysis (ASCA) a multivariate extension of ANOVA analysis was performed to allow interpretation of the variation induced by the different factors including age, sex, case, and their interaction. Subsequently, supervised partial least square discriminant analysis (PLS-DA) was performed to discriminate between the samples groups (PT1D/P1Ab/CTR) using combination of various categorical variables (e.g. PT1D vs. CTR, PT1D vs. P1Ab) in a specific age cohort. The area under the receiver operating characteristic (AUROC) were calculated to evaluate the ability of discrimination model to correctly classify the samples. All PLS-DA models were cross-validated by splitting the dataset into test sample subsets for each sub-validation experiment using random subset with ten data splits and twenty repeats. Furthermore, supervised multilevel partial least square discriminant analysis (ML-PLSDA), a multivariate extension of paired t-test was performed to discriminate between matched paired of samples, for instance the before and after seroconversion samples. ML-PLSDA models were validated using double cross model validation (CMV) 48 with twenty repeats. The mat lab code for ML-PLSDA is open access and available at http://www. bdagroup.nl/.
For univariate analysis, lipid concentrations were compared using both parametric and non-parametric test as per requirement. Wilcoxon rank-sum test was performed for comparing the two study groups of samples (e.g. PT1D vs. P1Ab) in a specific age cohort assuming data violates heteroscedasticity. For comparison one sample per subject, closest to the age within the time window, has been used in each test. Paired t-test was performed for the matched groups of samples (e.g. before vs. after seroconversion). The resulting nominal p-values were corrected for multiple comparisons using Storey approach 49 if the number of tests in a study were large (>10) or Benjamin and Hochberg approach 50 if the number of tests were small (<10). The adjusted p-values < 0.1 (q-values) were considered significantly different among the group of hypotheses tested in a specific age cohort. All of the univariate statistical analyses were computed in MATLAB 2017b using statistical toolbox. The fold difference was calculated by dividing the mean concentration of a lipid species in one group by another, for instance mean concentration in the PT1D by the mean concentration in P1Ab, and then illustrated by heat maps. The individual lipids levels were visualized as box plot using GraphPad Prism 7 (GraphPad Software, La Jolla, CA, USA). Data availability. The lipidomics data described in this study are deposited at the MetaboLights database 51 with the acquisition number (MTBLS620). The associated meta-data were captured using the ISA-creator package available from the MetaboLights. All the data supporting the findings of this study are available from MetaboLights database or from the corresponding authors on reasonable request.
Ethical approval and informed consent. The ethics and research committee of the participating university and hospital at University of Tampere, Tampere Finland, approved the study protocol. The study was conducted according to the guidelines in the Declaration of Helsinki. Parent and/or legal guardian provided written informed consent for participation in the study.