The study was approved by the Institutional Bioethical Committee (RNN/71/14/KE) and patient recruitment was performed between 30.10.2014 and 07.10.2016 in the Department of Pediatrics, Diabetology, Endocrinology and Nephrology of the Medical University of Lodz. The study protocol and its scientific purposes were explained to the patients and their legal guardians. Patients meeting the inclusion criteria were asked to participate in the study. Written consent to participate in the study was obtained from the patients and/or their legal guardians. We prospectively recruited patients who met the following inclusion criteria: age between 2 and 19 years old; diagnosis of type 1 diabetes according to ISPAD 2014 criteria (Rubio-Cabezas et al., 2014); lack of known concomitant genetic disorders; and diabetes treatment with insulin alone. Additionally, group-specific inclusion criteria were as follows: for the HG group—glucose concentration below 70 mg/dl during the episode and being unconscious or seizures; for DKA group—pH < 7.30 and HCO₃ ⁻ <15 mmol/L at diabetes diagnosis; for NDM group—pH > 7.35 and HCO₃ ⁻ >21 mmol/L at diabetes diagnosis; and for EDM group—diabetes duration above 6 months and lack of DKA or HG episodes within last 3 months. Clinical data were collected from the hospital information system. All patients within the batches were age- and sex-matched.

The sample collection protocol is shown in Figure 1A. For the HG group, the first sample was collected right after hospital admission (0 h), the second sample was collected 12 h after the admission, and the third sample was collected 48 h after the admission. For the DKA and NDM groups, the first samples were collected at hospital admission (0 h), the second—after 24 h admission, and the third ones—after 72 h after admission. Different protocols of the sample collection for HG and DKA/NDM patients arose from different dynamics of metabolic disturbances typical for the evaluated groups and the fact that HG patients were admitted after the episode and DKA/NDM were admitted with an ongoing disturbance and different hospitalization time: 2–3 days for the HG group and up to 2 weeks for the DKA–NDM group. Also, a limiting factor was the availability of an intravenous line in patients which is usually kept up to 3 days. For the EDM group only one sample was collected during routinely performed control laboratory tests as we did not assume a dynamic change in their metabolomic profile. Due to ethical issues, blood samples were collected only from already installed intravenous catheter or during routinely performed blood collections. Each time, a sample of 4 ml was collected from the patient to the vial with a clot activator (silica particles) (Becton Dickinson, New Jersey, NJ, United States). In order to standardize the sample processing protocol, we used vials from the same manufacturer (series 8516328 BD/REF 369032/2015–08/4113194). Between 60 and 120 min after the blood collection, serum isolation was performed by centrifugation for 10 min with 800 g. After serum separation, the biological material was frozen in -80°C degree until metabolic fingerprinting.

On the day of the metabolomic analysis, samples were thawed on ice. Protein precipitation and metabolite extraction were performed by vortex-mixing (for 1 min) one volume of the serum sample with four volumes of freeze cold (–20°C) methanol/ethanol (1:1) mixture containing 1 ppm of zomepirac (IS). After extraction, the samples were stored on ice for 10 min and centrifuged at 21,000 × g for 20 min at 4°C. The supernatant was filtered through a 0.22-µm nylon filter into glass vials.

Quality control (QC) samples were prepared by mixing an equal volume of all samples. The obtained mixture was prepared following the same procedure as the rest of the samples. QC samples were used at the beginning of the analysis (10 injections), in order to obtain the system stability, and then after every nine samples, in order to control the system stability and productiveness of the sample measurements.

Blank samples included mixtures of methanol/ethanol (1:1) that were prepared in the same way as biological samples. Zomepirac sodium salt (used as the internal standard (IS), arachidonic acid, docosahexaenoic acid, lysophosphatidylethanolamine (LPE) 18:0, phosphatidylethanolamine (PE) 16:0/22:6, LS-MS-grade acetonitrile, methanol, LC-grade formic acid, and ethanol were purchased from Sigma-Aldrich Chemie GmbH (Steinheim, Germany). The API-TOF reference mass solution kit (G1969-850001) and tuning solutions, ESI-L low-concentration tuning mix (G1969-85000), and ESI-TOF Biopolymer Analysis reference masses (G1969-850003) were purchased from Agilent Technologies (Santa Clara, CA, United States). Purified water was obtained using the Milli-Q Integral three system (Millipore SAS, Molsheim, France) in order to prepare solution A.

Metabolic fingerprinting was performed using liquid chromatography coupled to the tandem (quadrupole and time of flight, Q-TOF) mass spectrometry (LC-MS) system following the previously described protocol (Daniluk et al., 2019). Samples from the HG and DKA groups were analyzed in the separate batches to ensure data quality.

Analyses were performed in compliance with the current standards in force for the metabolomic measurements based on mass spectrometry and quality control methodology (Dunn et al., 2011; Godzien et al., 2015a).

Samples were randomly analyzed using the LC-MS system consisting of 1290 Infinity LC with a degasser, two binary pumps, and a thermostated (4°C) autosampler coupled to a 6550 iFunnel Q-TOF-MS detector (both Agilent Technologies, Santa Clara, CA, United States). Analyses were performed in positive (ESI+) and negative (ESI-) ion modes, whereby 1 µL of the sample was injected into a thermostated (60°C) Zorbax Extend- C18 RRHT (2.1 × 50 mm, 1.8-µm particle size, Agilent Technologies) chromatographic column. The flow rate was 0.6 ml/min with solvent A (water with 0.1% formic acid) and solvent B (acetonitrile with 0.1% formic acid). The chromatographic gradient started at 5% of phase B for the first minute, followed by an increase of phase B to 80% (from 1 to 7 min) and to 100% (from 7 to 11.5 min). After reaching 100%, the gradient returned to their initial conditions (5% phase B) in 0.5 min, which were kept from 12 to 15 min.

The mass spectrometer operated in the full scan mode from 50 to 1,000 mass (m/z). The capillary voltage was set to 3,000 V for the positive and 4,000 V for the negative ionization mode. Nozzle voltage was 1,000 V. The drying gas flow rate was 12 L/min at 250°C and gas nebulizer at 52 psig; the fragmentor voltage was 250 V for the positive and negative ionization modes. Data were collected in the centroid mode at a scan rate of 1.5 scan per second. Accurate mass measurements were obtained by means of a calibrant solution delivery using a dual-nebulizer ESI source. A calibrating solution containing reference masses at m/z 121.0509 (protonated purine) and m/z 922.0098 (protonated hexakis [(1H,1H,3H-tetrafluoropropoxy) phosphazine or HP-921] in the positive ion mode or m/z 119.0363 (proton abstracted purine) and m/z 966.0007 (formate adduct of HP-921) in the negative ion mode was continuously introduced by an isocratic pump (Agilent, Santa Clara, CA, United States) at a flow rate of 0.5 ml/min (1:100 split).

The raw data collected by the analytical instrumentation were cleaned of background noise and unrelated ions by the molecular feature extraction (MFE) tool in the Mass Hunter Qualitative Analysis Software (B.07.00, Agilent, Santa Clara, CA, United States). The MFE creates a list of all possible components described by mass, retention time (RT), and abundance.

The limit for the background noise for data extraction by MFE was set to 1,500 and 1,000 counts for the positive and negative ion modes, respectively. To identify co-eluting adducts of the same feature, the following adduct settings were applied: + H, + Na, and + K in the positive ion mode and −H, + HCOO, and + Cl for the negative ion mode. Dehydration neutral losses were also allowed in both ionization modes. Sample alignment and data filtering were performed using Mass Profiler Professional 12.6.1 (Agilent, Santa Clara, CA, United States). Parameters applied for the alignment were 1% for RT and 15 ppm for the mass variation. In the quality assurance (QA) procedure, metabolic features detected in >80% in QC samples with the coefficient of variation (CV) < 20% were kept for further data treatment.

Accurate masses of features were searched against the METLIN (http://metlin.scripps.edu), KEGG (https://www.genome.jp/kegg/pathway.html), LIPIDMAPS (http://lipidmaps.org), and HMDB databases (http://hmdb.ca), which were simultaneously accessed by a CEU Mass Mediator (http://ceumass.jpg.uspceu.es/mediator/). The identity of the metabolites was confirmed by matching the experimental MS/MS spectra to the MS/MS spectra from databases or fragmentation spectra and retention time obtained for the metabolite’s standard. Experiments were repeated with identical chromatographic conditions to the primary analysis. Ions were targeted for collision-induced dissociation (CID) fragmentation on the fly based on previously determined accurate mass and retention time. Phospholipids were identified based on a previously described characteristic fragmentation pattern (Godzien et al., 2015b). Details on metabolite identification are provided in Supplementary Table S1. Internal standard verification of two metabolites for which IS was commercially available is shown in Supplementary Table S2.

A quality assurance procedure was performed and the metabolomic features with RSD (relative standard deviation) for QC samples ≥20% or detectable in more than 1/3 blank samples were filtered out. Due to the high variance of a total signal obtained from each sample (probably due to samples dilutions caused by ongoing fluid therapy of the patients), the levels of themetabolomic features were divided by the sample total signal and multiplied by 100,000,000. This should make the metabolomic feature levels more comparable between samples with different dilutions. One sample from the NDM group (X62-5—second time point) was found to be an outlier in the negative ionization mode in the HG batch due to very low total signal, and most metabolic features undetectable—probably due to a problem with injection to the analytical system (Supplementary Figure. S1). It was excluded from further analysis.

For the statistical analysis, metabolomic features occurring in at least 80% of samples from each group were selected. A univariate comparison was performed with ANOVA and Benjamini–Hochberg for multiple comparison correction. For the DKA series, 48h and 72 h samples were compared with the NDM group samples from respective time points. DKA samples from the first time point were omitted here due to an ongoing DKA state that greatly affected the results of this analysis. For HG series, all time point samples were used to compare them with the EDM group. In this analysis, each time point was kept as a separate group decreasing statistical power thus for the post-hoc analysis, features with FDR<0.15 were selected. A Bonferroni-adjusted t-test was used for post-hoc pairwise comparisons. In the univariate analysis, for the DKA batch, the NDM group was used as control, while for the HG batch, the EDM group was used as control.

An advanced statistical analysis was performed with the Metaboanalyst 4.0 tool (Chong et al., 2018) with data filtered in the same way as for a simple statistical analysis. Additionally, missing data were filled with the K-NN (nearest neighbor) method and were scaled with the Pareto algorithm. An MS peak to Pathway analysis was performed with the GSEA algorithm and with the MFN database. Two OPLS-DA (orthogonal partial least square discriminant analysis) models were built for each study group. The first one was to discriminate between study groups from EDM samples, and the second one to discriminate from all NDM samples. Metabolomic features meeting the following criteria: p[1]>|0.4| and p(corr)[1]>|0.2| from both models were selected.

For features meeting the selection criteria in the simple and advanced analyses, an identification attempt was made.

Collected clinical data were compared between the groups by parametric or non-parametric tests based on data distribution and homoscedasticity. Metabolite levels were log-transformed with base 10. Correlations between the metabolite levels and other clinical data were performed with Pearson or Spearman correlation—depending on the variable distribution. Diagnostic utility of identified metabolites was performed with the ROC curve analysis. Sensitivity and specificity along with 95%CI (confidence interval) were calculated with VassarStats (http://vassarstats.net/). A clinical data analysis was carried out with STASTISTICA 13.1 (TIBCO Software, Palo Alto, CA, United States).

For the DKA analysis, a batch of samples was collected from: 20 patients hospitalized due to DKA (with three samples collected per patient 0h-24h–72h, N_samples = 60), 10 patients with newly diagnosed type 1 diabetes mellitus (NDM, with three samples collected per patient 0h-24h–72h, N_samples = 30), and 10 patients with established diabetes mellitus (EDM, one sample collected per patient N_samples = 10). For the HG batch, we recruited 10 patients hospitalized after the episode of HG (with three samples collected per patient 0h-12h–48h, N_samples = 30), 25 patients with EDM (1 sample collected per patient N_samples = 25), and 15 patients with NDM (3 samples collected per patient 0h-24h–72h, N_samples = 45). Data collection protocol for each group is shown in Figure 1A. Study group characteristics for the DKA batch is shown in Table 1 and for the HG batch in Table 2.

The protocol of data filtering in both ionization modes for the DKA batch and the HG batch is shown in Supplementary Figure S2. After technical filtering in the DKA batch, 543 metabolomic features (248 and 295 in positive and negative, respectively) remained, while in the HG batch, 733 features were deemed eligible for further analysis (359 and 374, respectively) (Figure 1B,C).

Initially, using all the metabolomic features passing technical data processing, we searched for pathways disturbed by DKA/HG beyond the acute phase of the episodes (Supplementary Tables S3, S4). Thus, we performed the MS peak to Pathway analysis by using the Metaboanalyst tool. By those means, we found that after the DKA episode following pathways were enriched: “xenobiotic metabolism”, “anti-inflammatory metabolites formation from EPA”, “arachidonic acid metabolism”, and “prostaglandin metabolism” and “glycerophospholipid metabolism” were diminished (Figure 1D). After the HG episode, we detected only down-regulated pathways and they were associated with “hormone biosynthesis”, “hexose phosphorylation”, “sphingolipid metabolism”’, and “arachidonic acid metabolism” (Figure 1E).

For the DKA batch, after overlapping results of univariate and advanced statistical modeling (both OPLS-DA models), 22 metabolomic features were selected as the fingerprint of the DKA episodes (Figure 2A). Seven features were successfully identified as four metabolites—two sphingomyelins with a higher level after the DKA episode versus comparative groups: SM (34:0) (Figure 2B) and SM (d18:0/15:0) (Figure 2C) as well as LPC (18:1) (2 features) (Figure 2D) which had a lower level during the episode but higher in the post-episode period. One metabolite had a lower level after DKA—LPC (18:0) (2 features) (Figure 2E). A group-specific time point-paired depiction of the identified metabolite level is shown in Supplementary Figure S3.

For the HG batch, after overlapping results of univariate and advanced statistical modeling (both OPLS-DA models), nine metabolomic features were selected as a fingerprint of the HG episodes (Figure 3A). From those, four features were identified as four metabolites—3 with higher after the HG episode versus comparative groups: two lysophosphatidylethanolamine: (LPE) (20:3) (Figure 3B) and LPE (18:2) (Figure 3C) as well as LPC (18:2) (Figure 3D); and one with lower level in HG group - oxy-phosphatidylocholine [PC O-(34:3)] (Figure 3E). A group-specific time point-paired depiction of the identified metabolite level is shown in Supplementary Figure S4.

In the DKA batch, two sphingomyelins strongly correlated with one another (r = 0.94, p < 0.0001) (Supplementary Figure S5A) and their correlation with other metabolites: strong positive correlation with LPC (18:0) and a lack of correlation with LPC (18:1) were almost identical (Supplementary Figures S5C–F). Thus, for further analysis, we showed only results for one of them - SM (34:0). SMs’ strong correlations may indicate a similar mechanism leading to their increased level after the DKA episode, however two LPCs: (18:0) and (18:1) poorly correlated with each other (r = 0.18, p = 0.0180, Supplementary Figure S5B) suggesting different origins of their serum levels.

In the HG batch, LPC (18:2) and LPE (18:2) were strongly, positively correlated between each other (r = 0.80, p < 0.0001, Supplementary Figure S6A), as well as LPE (18:2) and LPE (20:3) (r = 0.49, p < 0.0001, Supplementary Figure S6B). Additionally, LPE (20:3) and LPC (18:2) positive, weakly correlated with one another (r = 0.37; p = 0.0003, Supplementary Figure S6C). Negative correlations were detected between: PC O-(34:3) and LPE (20:3) (R = -0.36; p = 0.0003, Supplementary Figure S6D), PC O-(34:3) and LPE (18:2) (R = -0.24, p = 0.0157, Supplementary Figure S6E), and between LPC (18:2) and PC O-(34:3) which was at the borderline of statistical significance (R = -0.17; p = 0.0902, Supplementary Figure S6F).

Clinical data were clustered to show common changing patterns and help the interpreting correlations with metabolite levels (Figure 4).

For the DKA markers, during the first 2 days, we found that sphingomyelins correlated strongly and negatively with all four DKA markers measured at admission (pH, pCO₂, HCO₃ ⁻ and BE) (Figure 4A). This indicates that sphingomyelin changes were associated with the depth of ketoacidosis and thus strengthens their role as DKA biomarkers. Additionally, we observed a positive correlation between LPC (18:0) and MPV, which was previously shown to be an indicator of the dehydration level during DKA (Małachowska et al., 2020). Also, a negative correlation between both sphingomyelins and HDL were found statistically significant.

For the HG markers, we noticed that positive markers correlated positively with hematocrit, hemoglobin, and RBC, and the negative marker correlated negatively with them (Figure 4B). Those correlations strengthened over time with the weakest being observed in time point just after the episode. This might indicate that red blood cells might be associated with the distorted metabolism of those metabolites in the days following the episode of HG.

To check if the selected metabolites could be used as biomarkers of a past episode of the DKA/HG, ROC curve analysis was performed (Supplementary Figures S7–8). The best biomarker of the past DKA episode (after normalization of the blood–gas test, so at least after 24 h since admission) would be LPC (18:1) as it exceeded 0.8 threshold for AUC for the 48 h time point and almost exceeded it for the 72 h time point (AUC = 0.789). For HG, the best diagnostic criteria in all time points were attributed to LPE (20:3) with all AUC values exceeding 0.8 and reaching 100% sensitivity.