The Footprint of Type 1 Diabetes on Red Blood Cells: A Metabolomic and Lipidomic Study

The prevalence of diabetes type 1 (T1D) in the world populations is continuously growing. Although treatment methods are improving, the diagnostic is still symptom-based and sometimes far after onset of the disease. In this context, the aim of the study was the search of new biomarkers of the disease in red blood cells (RBCs), until now unexplored. The metabolomic and the lipidomic profile of RBCs from T1D patients and matched healthy controls was determined by NMR spectroscopy, and different multivariate discrimination models were built to select the metabolites and lipids that change most significantly. Relevant metabolites were further confirmed by univariate statistical analysis. Robust separation in the metabolomic and lipidomic profiles of RBCs from patients and controls was confirmed by orthogonal projection on latent structure discriminant analysis (OPLS-DA), random forest analysis, and significance analysis of metabolites (SAM). The main changes were detected in the levels of amino acids, organic acids, creatine and phosphocreatine, lipid change length, and choline derivatives, demonstrating changes in glycolysis, BCAA metabolism, and phospholipid metabolism. Our study proves that robust differences exist in the metabolic and lipidomic profile of RBCs from T1D patients, in comparison with matched healthy individuals. Some changes were similar to alterations found already in RBCs of T2D patients, but others seemed to be specific for type 1 diabetes. Thus, many of the metabolic differences found could be biomarker candidates for an earlier diagnosis or monitoring of patients with T1D.


Introduction
Diabetes mellitus is the first non-communicable disease that the United Nations recognized as a 21st-century pandemic problem [1]. Diabetes is a group of diseases mainly characterized by high blood glucose levels. Type 1 diabetes (T1D) and type 2 diabetes (T2D) both involve hyperglycemia, but are disorders with different patterns, mechanisms, clinical presentations, and disease progressions. Nevertheless, once hyperglycemia occurs, patients with all forms of diabetes can develop similar chronic complications, in different rates of progression and form [2].
T2D accounts for about 80-90% of the cases of diabetes and is due to a progressive loss of β-cell mass and insulin secretion, on the background of insulin resistance (IR), low-grade inflammation, and oxidative stress [3]. In 2019, 463 million adults (20-79 years) were living groups as follows: 20 patients with T1D matched by age, sex, and body mass index (BMI) with 20 healthy individuals. All patients signed the informed consent form before their inclusion in the study.
The inclusion criteria were the following: (a) age > 18 years, (b) T1D, (c) written informed consent form. The exclusion criteria were the following: (a) any other type of diabetes, except for T1D (T2D, MODY, etc.), (b) tobacco use in the last 12 months, (c) alcohol intake of more than 30 g in males and 20 g in females, (d) hemoglobinopathies, (e) hematologic disease, (f) any other disease or treatment that might influence the RBCs' metabolism (chronic kidney disease, cardiac valve prosthesis, chemotherapy, anticoagulants, etc.). All diabetic patients had positive islet autoantibodies of anti-GAD-65 (>50 IU/mL).
This was a pilot study. At present, there are no available data regarding the impact of T1D on the RBCs' metabolism; therefore, we were not able to perform a sample size calculation. Thus, the size has been selected by taking into account ethical considerations and the pilot studies in T2D [25,26].
All the subjects included in the study underwent a complete medical history review, anthropometric evaluation, and biochemical analysis. Anthropometrical parameters of weight (kg), waist (cm), and height (m) were measured just before the collecting of the blood samples. Biochemical analysis was performed using standardized and validated routine methodologies on the blood samples collected at fasting conditions at the biochemistry core facilities of Vall d'Hebron University Hospital. Plasmatic concentrations of glucose (mg/dL), glycosylated hemoglobin (HbA1c) (%DCCT), high-density lipoprotein (HDL) (mg/dL), low-density lipoprotein (LDL) (mg/dL), triglycerides (TG) (mg/dL), anti-GAD65 (IU/mL), and insulin (mU/mL) were measured in all the subjects. Body mass index (BMI) kg/m 2 was calculated using the following formula: body weight (kg)/height 2 (m) [27].

Peripheral Blood Samples
Blood samples were collected followed a similar protocol as in previous studies [25,26]. Briefly, peripheral blood was collected under fasting conditions, stored at 4 • C, and processed within the first hour. Then, 5 mL of peripheral blood were poured into a quartz tube with 10 mL of Ficoll until 2 phases were separated by gravity, around 30 min. Then, the pellet was washed twice with 10 mL of cool Hank's Balanced Salt Solution in a centrifuge at 200× g and 4 • C for 10 min without breaks. Cell counting was performed on a Neubauer chamber and the purity was tested with flow cytometry. For storage, 200 µL of ice-cold methanol were added per 10 million cells for quenching and the samples were frozen directly at −80 • C.

Extraction of Metabolites and Lipids
Metabolites were extracted from RBCs as described previously [25]. Briefly, frozen samples were allowed to thaw for 5 min on ice. Then, 800 µL of chloroform at 4 • C were added to each sample tube. After 10 min, the samples were homogenized, resuspended with a pipette, and transferred to a bigger plastic tube. Then, samples were placed in liquid nitrogen for 1 min and then allowed to thaw on ice for 2 min, for uniform cell membrane breakage. This step was performed two more times. Afterwards, 1250 µL of distilled water and 1250 µL of chloroform were added and the sample vortexed. Samples were then centrifuged at 13,000× g for 20 min at 4 • C to separate the phases. The upper phase contained polar metabolites in a mixture of water/methanol and was lyophilized overnight. The lower chloroform phase containing the lipids was evaporated under nitrogen flux. Extracts were stored at −80 • C until NMR sample preparation and analysis. For NMR sample preparation, the frozen samples were allowed to thaw for 5 min on ice. Then, 550 µL of NMR buffer (100 mM Na 2 HPO 4 in D 2 O at pH 7.4) containing 0.1 mM deuterated trimethylsilylpropanoic acid as an internal standard were added to the polar samples and transferred to a 5 mm NMR tube. For the preparation of the lipidomics samples, 600 µL of deuterated chloroform containing 0.003% v/v tetramethylsilane was added to the lipidic samples and transferred to a 5 mm NMR tube. All samples were analyzed the same day and stored at 4 • C until the analysis.

NMR Experiments
NMR spectra were recorded as described previously [25] at 27 • C on a Bruker AVII-600 spectrometer using a 5 mm triple resonance cryoprobe and processed using TopSpin version 3.2 software (Bruker BioSpin). Before starting sample measurement, the temperature was calibrated at 27 • C with a 99.8% MeOD sample. Further, the parameters of the spectrometer were optimized to ensure an optimal resolution and sensitivity (width at half-height ≤ 1 Hz), with a standard sample containing 2 mM of sucrose, 0.5 mM of DSS (sodium trimethylsilylpropanesulfonate), 2 mM of NaN 3 in 90% H 2 O, and 10% D2O. 1 H 1D NOESY NMR spectra were acquired with 256 free induction decays (FIDs). A 4s relaxation delay was incorporated between FIDs and the water signal was eliminated with presaturation in the aqueous samples.
The FID values were multiplied by an exponential function with a 0.5 Hz line broadening factor for an optimal baseline correction, with 64K data points digitalized over a spectral width of 30 ppm. Quality control samples, consisting of a mixture of lactate, creatine, citric acid, phenylalanine, and glucose, were run in regular intervals to confirm the reproducibility of the measurements.

Data Analysis
1 H-NMR signals from the spectra were assigned with the help of previous data [25,28] and the spectral databases Human Metabolome Database and Biological Magnetic Resonance Bank [29,30]. In ambiguous cases, the assignment was confirmed by spiking the spectra with reference compounds. For metabolite quantification, spectra were normalized to total intensity. This allowed us to work with normalized concentration values and avoid the detection of changes related exclusively to differences in total sample amount and experimental error during the extraction process. The predominant glucose signals and the solvent signals were excluded from normalization. Optimal integration regions were defined for each metabolite, an attempt being made to select the signals without overlapping (Table S1). Integration was performed with global spectrum deconvolution in MestreNova 12 (Mestrelab Research S. L.). The datasets are available in the Zenodo repository [31].
For modelling, metabolite integration data were Pareto-scaled (each value being divided by the square root of the standard deviation of each variable) and mean-centered, for an easier interpretation of the data and to also take into account the variation of small signals. Orthogonal projection on latent structure discriminant analysis (OPLS-DA) and projection on latent structure (PLS) regression were performed with SIMCA-P 14.0 (Umetrics, Sweden). OPLS-DA and PLS models were described with R2Y(cum) (representing the cumulative SS of all the y-variables explained by the extracted components) and Q2(cum) (describing the cumulative Q2 for all y-variables for the extracted components). Score plots (representing scores related to prediction Y t(1) versus scores related to the first orthogonal component t0) and S-plots (plotting the modelled covariation (p(1)) versus the modelled correlation (p(corr)) [32] were used to detect significant metabolites and lipids in the OPLS-DA models. The glucose signal, due to its high intensity and correlation to plasmatic glucose, was excluded from OPLS-DA analysis. OPLS-DA and PLS models were validated in SIMCA-P by a 100-times permutation (where the overfit of the model is measured by the intercept of the regression line of the correlation coefficient between the original y-variable and the permuted y-variable on the x-axis versus the cumulative R2 and Q2 on the y-axis) and analysis of variance testing of cross-validated predictive residuals (CV-ANOVA) [33]. CV-ANOVA is a significance test for the Q2YCV cross-validation using the F-distribution, based on an ANOVA assessment of the cross-validatory (CV) predictive residuals of the model. In PLS analysis, metabolites that were relevant for the correlation were identified by VIP (variable-importance-projection) values >1.
Significance analysis of metabolites (SAM) analysis was performed with Metaboanalyst [34] applying the siggenes package [35]. To each variable, a significance score was assigned based on its change relative to the standard deviation of repeated measurements. For a variable with scores greater than an adjustable threshold, its relative difference is compared to the distribution estimated by random permutations of the class labels. For each threshold, a certain proportion of the variables in the permutation set will be found to be significant by chance. The proportion is used to calculate the FDR. In our analysis, a Delta value of 1.4 was defined to control FDR.
Random Forest (RF) analysis was performed with Metaboanalyst [34] based on the randonForest package [36]. RF uses groups of classification trees, each of which is grown by random feature selection from a bootstrap sample at each branch. Class prediction is based on the majority vote of the ensemble. For initial tree construction, one-third of the samples are left out for validation. The OOB (out-of-bag) data are then used as a test sample to obtain an unbiased estimate of the classification error (OOB error). Variable importance is calculated by the measure of the increase of the OOB error when it is permuted.
Hierarchical clustering was performed with Metaboanalyst with the hclust function in package stat and the result represented as a heatmap.
Box plots and t-tests were performed with normalized concentration data in Graphpad Prism 8. Only metabolites that showed high scores in all three models (OPLS-DA, Random Forest, and SAM) were considered for univariate analysis.
For the biological interpretation of the results and the identification of metabolic pathways, the Kegg Data Base and MetPA (Metaboanalyst) were used. Table 1 summarizes the main characteristics of the subjects and controls included in the study. A metabolomic and a lipidomic analysis of RBCs of all individuals was performed by proton nuclear magnetic resonance ( 1 H-NMR). As a result, 68 different hydrophilic metabolites and 12 different types of lipid groups were identified and quantified (Supporting Tables S1 and S2). Three different discrimination methods were applied in order to analyze significant differences between both groups: OPLS-DA (Orthogonal Projections to Latent Structures Discriminant Analysis), Random Forest Classification, and SAM (Significance Analysis of Metabolites). Significant discrimination between healthy controls (CT) and T1D patients could be obtained for the three model types of the metabolomic and lipidomic profiles (Figure 1).

Results
In addition, a heatmap analysis was performed to have an overview of the main metabolic changes in the patient and control groups. As can be seen in Figure 2A, a similar tendency in metabolite levels in T1D patients is detected (yellow color), while controls have a higher dispersion (blue color). As a general tendency, T1D patients seem to have lower levels of amino acids and organic acids.
Concerning lipid metabolism, T1D patients tended to have lower levels of unsaturated fatty acids ( Figure 2B). This tendency was, however, less pronounced for polyunsaturated fatty acids (PUFA). Additionally, phosphatidylcholines were reduced as well as phospholipids in general. On the other hand, a decrease of lipid chain length could be detected in  In addition, a heatmap analysis was performed to have an overview of the main metabolic changes in the patient and control groups. As can be seen in Figure 2A, a similar tendency in metabolite levels in T1D patients is detected (yellow color), while controls have a higher dispersion (blue color). As a general tendency, T1D patients seem to have lower levels of amino acids and organic acids. In order to make a robust metabolite selection, metabolites that had high scores in all three model types (Figure 1) were further submitted to univariate analysis. The resulting significant metabolites are represented in Figure 3, and can be divided into different metabolite groups. As can be seen, organic acids and metabolites related to lipid and choline metabolism seem to be the most affected. detected in T1D patients, determined by the ratio of lipid methyl groups (Lipid CH3-, increase) versus linear lipid chains (Lipid-CH2, decrease).
In order to make a robust metabolite selection, metabolites that had high scores in all three model types (Figure 1) were further submitted to univariate analysis. The resulting significant metabolites are represented in Figure 3, and can be divided into different metabolite groups. As can be seen, organic acids and metabolites related to lipid and choline metabolism seem to be the most affected. In order to detect which of the metabolic changes were directly correlated with plasma glucose levels in fasting conditions, we performed PLS regression analyses of the metabolomic and the lipidomic profile versus glucose concentrations. As a result, we obtained a good correlation of glucose with several metabolites ( Figure S1A) and a moderate correlation with some blood lipids ( Figure S1B). The metabolites and the lipids that were relevant for the correlation (VIP > 1) are listed in Tables S3 and S4. Most metabolites that were significantly altered in RBCs of T1D patients also had a correlation with plasma In order to detect which of the metabolic changes were directly correlated with plasma glucose levels in fasting conditions, we performed PLS regression analyses of the metabolomic and the lipidomic profile versus glucose concentrations. As a result, we obtained a good correlation of glucose with several metabolites ( Figure S1A) and a moderate correlation with some blood lipids ( Figure S1B). The metabolites and the lipids that were relevant for the correlation (VIP > 1) are listed in Tables S3 and S4. Most metabolites that were significantly altered in RBCs of T1D patients also had a correlation with plasma glucose levels. Only glycerophosphocholine (GPC), inosinic acid (IMP), and creatine were not relevantly correlated. Regarding lipidomics, although the correlation was much weaker, all relevant lipid groups (lipid CH 3 , lipid CH 2 , and phosphatidylcholine) seemed to be related to plasma glucose concentration.
In order to correlate the metabolite changes with alterations of specific metabolic pathways, an enrichment pathway analysis was performed (Figure 4), showing that energy metabolism, fatty acid metabolism, and amino acid metabolism were the most affected metabolic pathways in T1D compared to healthy controls. not relevantly correlated. Regarding lipidomics, although the correlation was much weaker, all relevant lipid groups (lipid CH3, lipid CH2, and phosphatidylcholine) seemed to be related to plasma glucose concentration.
In order to correlate the metabolite changes with alterations of specific metabolic pathways, an enrichment pathway analysis was performed (Figure 4), showing that energy metabolism, fatty acid metabolism, and amino acid metabolism were the most affected metabolic pathways in T1D compared to healthy controls.

Discussion
In the present pilot study, we detected a clear difference in the metabolomic and lipidomic profile of patients with T1D and control subjects, proving that T1D has a specific

Discussion
In the present pilot study, we detected a clear difference in the metabolomic and lipidomic profile of patients with T1D and control subjects, proving that T1D has a specific impact on RBC metabolism and transport. Our group has previously reported that the metabolomics profile of RBCs is generally altered in T2D disease [25], mainly affecting redox and energy metabolism, as well as the transport of certain metabolites with physiological functions. As far as we know, this is the first study that shows a significant impact of T1D on the metabolomics and lipidomics profile of RBCs.
In our analysis, we detected a significant reduction of acetate in patients, as well as a tendency of lactate to decrease (Figure 2A). This could reflect a reduction in glycolysis, a simple pathway of glucose metabolism that is regulated by insulin secretion and sensibilization [37]. Previous studies reported increased plasmatic lactate and acetate levels in patients with diabetes [38]. Thus, the increment in lactate could be related to increased carbonic anhydrase activity in RBCs, transporting a higher amount of lactate into blood plasma [39]. Our group reported similar findings in patients with T2D [25], which could indicate that similar processes happen in T1D and T2D.
Other metabolites altered in RBCs were isoleucine and leucine, essential metabolites for the production and formation of hemoglobin and the production of RBCs [40]. This alteration is similar to that detected in T2D patients [25], where the BCCAs valine and leucine decreased. An opposite tendency was detected in plasma of diabetic patients where the concentration of BCAAs was increased in both T1D and T2D [41,42]. This inverse correlation could be related to the fact that RBCs with highly glycosylated hemoglobin, as occurs in both T1D and T2D diabetes, seem to have a lower BCAA uptake capacity. Interestingly, isoleucine was exclusively altered in RBCs of T1D patients, while significant changes in valine were only detected in T2D. The differences in the alteration of BCCAs depending of the type of diabetes could indicate a different pathophysiology associated with BCAAs in both types of diabetes, since this type of amino acids is transported by RBCs. Furthermore, valine and leucine could be RBCs' fingerprints of T2D and T1D, respectively. In addition, the alteration of the organic acids 2-hydroxyisovalerate and 3-mehtyl-2-oxovalerate, key compounds in BCAAs metabolism, were only decreased in T1D, indicating a higher affectation of BCAAs in this type of diabetes.
Previous studies reported changes in other kinds of amino acids associated with diabetes and pre-diabetes in many studies [43]. Interestingly, the changes in non-BCAA amino acids seemed to be less pronounced in RBCs of T1D patients than in T2D [25]. For instance, alanine, asparagine, glutamate, proline, threonine, and methionine seem not to be affected by T1D, while they change significantly in T2D. These altered amino acid levels may directly interfere with insulin-regulated glycogen synthesis, and the lower impact on amino acid concentration in RBCs may be related to a lower level of IR in T1D, reducing the enhancement of plasma concentration of these type of metabolites [44]. However, further studies are mandatory to clarify this assumption.
Another noteworthy result from our study was the increase of creatine and phosphocreatine in RBCs of T1D patients. Plasmatic creatine has been recently identified as a potential biomarker for mitochondria dysfunction and for T2D risk in the PREVEND study from Post et al [45]. Our results indicate that this increase could be more generalized in the organism, taking place in cells, and may be related to alterations in muscle or energy metabolism where creatine plays an important role. This is further supported by the increase of inosine monophosphate (IMP) showed in RBCs of T1D patients, a nucleotide related to energy and generated in muscles during exercise [46]. Interestingly, neither creatine nor IMP showed a direct correlation with blood glucose. Thus, further studies could prove their usefulness as hyperglycemia-independent biomarkers for T1D in early stages or non-fasting conditions. Interestingly, we did not detect significant changes in the levels of reduced and oxidized glutathione (GSH and GSSG), although we detected a tendency of GSH levels to decrease, and the glutathione pathway seemed to be affected in enrichment analysis ( Figure 4). Decreased GSH and increased GSSG levels have been previously detected in plasma [27,47,48] and RBCs [49,50] of T2D patients. This alteration was attributed to a decrease in the activities of the enzymes involved in GSH synthesis [49], a decrease in the transport rate of oxidized glutathione (GSSG) from RBCs [51], an enhanced sorbitol pathway [52], and the increase of oxidative stress produced by the elevated generation of ROS linked to T2D [53]. Our result was coherent with the fact that a decrease of GSH in blood of T1D patients was detected less frequently [47,54]. Thus, the glutathione pathway seems to be less affected by T1D.
We also did not detect significant differences in 2,3-Bisphosphoglyceric acid (2,3-BPG) levels, which have been described for T2D [25]. 2,3-BPG regulates oxygen absorption by hemoglobin, and its alteration can be directly associated to hypoxia and related to oxidative stress [55]. Thus, our study indicates lower hypoxia in organs and tissues related to RBCs and reactive oxygen species (ROS) levels in RBCs in T1D patients, in comparison with T2D. This result is in accordance with the lower impact on glutathione metabolism that we detected for T1D patients.
Regarding the lipidomic profile of RBCs, we found very significant differences between T1D patients and healthy controls. This result was in coherence with previous works that detected alterations in lipid metabolism associated with T1D and T2D [56][57][58]. Specifically, we detected a very clear decrease in lipid chain length in combination with an increase in lipid methyl groups, which could be an indication for lipid peroxidation, as was observed previously in peripheral blood mononuclear cells [59]. On the other hand, the increase in methyl lipid groups may also be related to the higher lipid accumulation of T1D patients due to the alterations of insulin release from pancreatic beta cells that inhibit the glucose uptake in cells [60]. In addition, we detected changes in the phospholipid phosphatidylcholine, as well as for the related metabolites phosphocholine and glycerophosphocholine (GPC), which are all important constituents of the cell membrane. As mentioned, diabetes negatively affects the membrane and the shape of RBCs [19,25,56], which is in coherence with these results. This alteration also was shown in RBCs of T2D patients [25].
Our study is only a first pilot study, and has several limitations that will have to be overcome in additional studies. For instance, although our study groups were sex-matched to compensate for metabolic differences between men and women, the low number of samples did not allow us to perform sex-specific studies to identify specific alterations for each sex group. Additionally, the effect of many other factors, such as age, comorbidities, or treatments, should be studied.
Nevertheless, we show robust differences in the metabolic and lipidomic profile of RBCs from T1D patients that provide a valuable starting point for the identification of biomarkers for an earlier diagnosis or monitoring of patients with T1D. Several of the detected metabolic changes are similar to alterations detected previously in the RBCs of T2D patients. Further studies are needed to confirm if these alterations are also related to the presence of IR in T1D patients. Informed Consent Statement: Informed consent was obtained from all subjects involved in the study.