Serum Metabolomic Signatures Indicate Oxidative Membrane Lipid Remodeling in β-Thalassemia

1. Introduction

β-Thalassemia is an autosomal-recessive hemoglobinopathy caused by deficient β-globin synthesis. Homozygous or compound heterozygous mutations produce transfusion-dependent (TD) disease with life-long anemia, whereas milder genotypes may be non-transfusion-dependent (NTD) with intermittent transfusion needs. Excess unpaired α-chains form hemichromes that catalyze reactive oxygen species, driving lipid/protein oxidation, membrane destabilization, apoptosis of erythroblasts, and ineffective erythropoiesis; peripheral hemolysis further aggravates anemia [1,2]. Beyond hematologic manifestations, chronic hemolysis, endothelial adhesion, inflammation, oxidative stress, and iron overload contribute to multisystem—especially cardiovascular—injury [3,4].

Metabolomics, by quantifying small-molecule intermediates and end products, offers a systems view of these processes [5,6]. Prior studies have reported altered organic and amino acids by NMR, targeted MS predictors of iron overload/oxidative damage, and GC-MS–defined serum pathway disturbances; hydroxyurea may partially normalize profiles [7,8,9,10,11]. Work in thalassemia models and patients highlights tissue-specific macrophage reprogramming, bile-acid perturbations, plasma/erythrocyte redox remodeling, and microbiome–metabolome links—including taxa that modify anemia and iron burden [12,13,14,15,16].

We therefore applied untargeted UHPLC–Orbitrap MS to serum samples from β-thalassemia patients versus healthy controls and explored associations between metabolites/pathways and clinical features (e.g., transfusion status, iron indices, MRI T2*, liver iron concentration, and therapy). Mapping these metabolic derangements aims to identify biomarkers and therapeutic targets that link the primary hemoglobin defect to multisystem disease. Such a transformation of hemoglobin-centered problem into a multisystemic disease will certainly help to discover new biomarkers and therapeutic strategies in β-thalassemia.

2. Materials and Methods

The study included 31 patients with β-thalassemia living on region of Northwestern Greece and followed at the pediatric and adult Thalassemia Units at the University Hospital of Ioannina. The diagnosis was based on the clinical, hematological, biosynthetic, and genetic studies. Medical files of the patients were retrieved, and detailed data were recorded and analyzed. Other parameters such as type of iron chelators, serum ferritin, cardiac MRI T2* and MRI liver iron concentration (LIC) were enrolled in the study. The demographic, laboratory and clinical data are shown in Table 1. Fourteen patients were male, and the age range was 12-67 years with a median 36 years. Eighteen had TD β-thalassemia and 13 NTD. Eight apparently healthy individuals matched for age and sex respectively, to patients, served as controls. All participants provided written informed consent, and the study was conducted in accordance with the Declaration of Helsinki.

Peripheral blood was drawn after an overnight fast and, for transfused patients, prior to the scheduled transfusion to minimize acute post-transfusion effects. Patients having any evidence of other chronic illnesses unrelated to thalassemia were excluded. Samples were allowed to clot, centrifuged to obtain serum, aliquoted to avoid repeated freeze–thaw, and stored at −80 °C until analysis. All steps were performed on ice or at 4 °C to limit ex vivo metabolic change and oxidation. Sterile pipettes and tips were used, and appropriate precautions were taken during the performance of assays (using a laminal flow instrument in a cell culture room) to protect samples from subsequent environmental contamination.

Acetonitrile, methanol, and formic acid were obtained from Fluka/Riedel-de Haën (Switzerland). Millipore’s Milli-Q Plus water purification system (Milford, MA, USA) was utilized to create high-quality purified water.

Serum samples were prepared following a specific protocol: initially, the samples were thawed on ice at 4 °C for 30 to 60 minutes. Then, 200 μL of each serum sample was placed into labeled 1.5 mL Eppendorf tubes, and 800 μL of methanol was added, followed by 15 seconds of vortexing. Afterward, the samples were centrifuged at 15800g for 15 minutes at room temperature to separate the protein precipitate. The resulting supernatant were transferred into two separate labeled 1.5 mL Eppendorf tubes. These were then evaporated in the GeneVac at 50 °C and stored (-80oC) until analysis. For HRMS analysis, 100 μL of 80:20 H2O:AcN was added to reconstitute the sample, vortexed for 15 seconds, centrifuged at 15800g for 15 minutes, and transferred to 200 μL inserts in appropriate screw-capped autosampler vials. Throughout these processes, the samples were kept on ice before each centrifugation, and all steps were conducted at a controlled temperature of 4 °C.

A pooled quality control (QC) sample was made to oversee the consistency of LC and MS performance across sample runs. This QC sample was formed by mixing equal small portions from all the samples in the experimental set before the evaporation process. Afterward, the combined mixture was evaporated until completely dry and managed consistently with the previously outlined procedures as a control throughout the experiments. Sample handling and analysis is depicted as a workflow in Figure 1.

Samples were analyzed using a Thermo Scientific ESI-LTQ-Orbitrap Discovery XL mass spectrometer connected to an Accela UHPLC system, also from Thermo Scientific. The UHPLC setup included an autosampler, vacuum degasser, binary pump, and a column with controlled temperature. Analysis was conducted using a Waters Corp. ACQUITY UPLC BEH C18 reversed-phase column. The system used a binary gradient solvent method with formic acid (0.1%v/v)-water (solvent A) and acetonitrile (solvent B). Sample analysis covered both positive (ESI+) and negative (ESI−) ion modes, maintaining a constant flow rate of 0.4 mL/min. The gradient method spanned 32 minutes with specific solvent ratios at different time points, while the column temperature remained at 40 °C and the autosampler tray temperature at 5 °C. The injection volume was 10 μL. In positive ion mode, the capillary temperature and voltage were adjusted to 320 °C and 40 V, respectively. The sheath gas flow was set at 40, aux gas flow at 8 arb. units, spray voltage at 3.30 kV, and tube lens voltage at 120 V. Conversely, for negative ion mode, the capillary temperature and voltage were set to 320 °C and −40 V, respectively. The sheath gas flow and aux gas flow were maintained at 40 and 8 arb. units, while the spray voltage was set to 2.70 kV and the tube lens voltage to −120 V. In both positive and negative ion modes, the analysis employed Fourier transform mass spectrometry (FTMS) in full scan ion mode, scanning within the range of mass-to-charge (m/z) ratio 100–1000 and a resolution of 30,000 FWHM, while spectra were acquired in centroid mode.

Before conducting the multivariate statistical analysis, the data underwent preprocessing using several tools. Xcalibur^® from Thermo Fisher Scientific, MZmine 2.10 available at mzmine.sourceforge.net. Microsoft Excel, and the web-based MetaboAnalyst suite [17] were utilized. Initially, Xcalibur^® converted the instrument’s native data files (.raw) to the cdf data format (.cdf) through the Xcalibur^® Xconvert program. Subsequently, MZmine 2.10 processed the data, employing baseline correction, peak detection, deconvolution, deisotoping, alignment, and gap filling procedures. The resulting peak list, displaying accurate mass - tR versus intensity, was exported as a .csv file to Microsoft Excel and further manipulated using commands such as CONCATENATE, ROUND, and TRANSPOSE for appropriate adjustments.

Next, the MetaboAnalyst workflow was used for statistical analysis. Within this analysis, several steps were taken: (a) missing values were addressed using the k-nearest neighbor (KNN) algorithm, (b) data filtering involved applying the interquartile range (IQR) to eliminate 25–40% of near-constant values throughout the experiment as a robust estimator of noise, and (c) row-wise normalization was conducted by averaging all samples (creating a pseudo-reference sample) in the control group. Additionally, either unit variance (UV) or Pareto scaling was explored to approximate a normal distribution as closely as possible.

Initially, Principal Component Analysis (PCA) was utilized to explore the inherent variations within the datasets, identify potential outliers through 2D scores plots, and assess the method’s quality by visually inspecting the QC clustering. Subsequently, Partial Least Square-Discriminant Analysis (PLS-DA) was employed in a sequential manner to unveil potential underlying classifications among the studied groups and to investigate the crucial features contributing to this clustering. The performance of the PLS-DA models underwent assessment using k-fold cross-validation (k = 10), leave-one-out cross-validation (LOOCV), and permutation testing (100 permutations) (18) R2 and Q2 values were employed for selecting the optimal parameters and evaluating the model’s predictive capacity. Variables with a Variable Importance in Projection (VIP) exceeding 1.5 were considered as highly influential in classifying the healthy controls versus patients groups. Additionally, orthogonal projection to latent structures–discriminant analysis (OPLS-DA) was employed using SIMCA P+ 10.5 (Umetrics, Umea, Sweden) and EZinfo 2.0 (Umetrics, Umea, Sweden). VIP scores and S-plots were constructed to assess the most influential metabolites in this scenario, while model quality was assessed based on R2X, R2Y, Q2, and permutation testing involving 100 random permutations.

Furthermore, an additional assessment of the generated outcomes involved conducting univariate data analysis through several methods: t-tests, fold change (FC) examination, and their amalgamation represented as a heatmap plot. Throughout these analyses, a consistent FC threshold of 2 was set, and the significance threshold was established at 95% (p-value = 0.05). To assess normality, the data underwent the Shapiro–Wilk, Kolmogorov–Smirnov, and Lilliefors tests using the ‘nortest’ R-package.

3. Results

3.1. Global Serum Metabolome and Data Structure

Untargeted UHPLC–Orbitrap MS detected 183 serum features across patients and controls. Based on combined multivariate and univariate analyses, 124 features were decreased, 54 were increased, and 5 were unchanged in patients versus controls

Two approaches were adopted for data scaling: UV (auto-scaling) and Pareto scaling for both negative and positive ion modes. UV scaling equalizes the importance of all metabolites but amplifies baseline noise. Conversely, Pareto scaling retains the original measurements to a greater extent but is responsive to significant changes [18]. As depicted in Figure 2 and Figure 3 for positive and negative ion modes, respectively, the samples showed close clustering in PCA score plots (data were not shown for UV scaling). However, the distinction between the two groups was not deemed satisfactory and therefore supervised methods, such as PLS-DA or OPLS-DA, were employed.

Both methods were able to distinguish between the two sets (controls versus patients), displaying distinctive clustering that suggested varying metabolic profiles in both modes both for UV and Pareto scaling. However, in these supervised models, that users categorize samples into pre-established groups, carries the risk of creating biased clustering or overfitting. Hence, permutation testing involving 100 permutations was used to evaluate the significance of the classification, utilizing the prediction accuracy during training as a statistical test. Τhe results of the Permutation test indicated that UV Scaling was over-fitted to the data since R2 and Q2 values derived from the permuted data exceeded the original values on the validation plot. This indicates that relying on the model using UV Scaling might not be dependable for subsequent statistical analysis, as it fails to precisely represent the inherent patterns within the data.

3.2. Feature Selection and Identification

The VIP-plot and S-plot in positive (Figure 4a,b) and negative (Figure 4c,d) mode were computed to pinpoint potentially crucial variables that contribute to the distinct clustering observed in the OPLS. VIP values (>1.0) were utilized to pinpoint the metabolites responsible for distinguishing between the two groups.

Additionally, to validate the outcomes related to the selection of VIPs from both modes (positive and negative), univariate data analysis was conducted using, t-test, FC monitoring, and their amalgamation in the form of a heatmap plot (Figure 5).

By merging the findings from both univariate and multivariate analyses, a compilation of features that differentiate the two groups was assembled for structural identification. The identification strategy involved comparing the selected features with online databases based on mass accuracy, isotopic patterns, RDBeq information, and potential low-energy or in-source fragmentation data for assessment. As a result, 20 metabolites (Table 2) were recognized via Orbitrap-based analysis.

Pathway analysis (MetaboAnalyst/KEGG) [17] showed significant enrichment for arachidonic acid metabolism (q = 0.022, impact = 0.433) and linoleic acid metabolism (q = 0.024, impact = 1.000) (Figure 6). Glycerophospholipid metabolism (q = 0.220, impact = 0.401) and porphyrin metabolism (q = 0.145, impact = 0.242) reached nominal significance but did not survive FDR correction. These patterns, together with LysoPC up-regulation at the feature level, support a PUFA (Polyunsaturated fatty acid)- and membrane-lipid–centered remodeling under oxidative stress, with eicosanoid precursor flux as a plausible downstream consequence. Primary bile-acid biosynthesis and α-linolenic acid metabolism were not significant after correction in this dataset.

3.4. Subgroup and Clinical Correlations

Within the β-thalassemia cohort, no significant differences were observed between TDT and NTDT groups in global multivariate structure or in individual metabolites. This lack of separation persisted after considering demographic, clinical, and laboratory covariates provided in Table 1.

4. Discussion

The present study utilized untargeted metabolomic analysis identified a set of significantly up-regulated metabolites that cluster into two biologically coherent groups: (A) phosphatidylcholine (PC)–related species along the glycerophospholipid axis and (B) conjugated bile acids belonging to primary/secondary bile-acid biosynthesis. Together, these findings support a unifying model in β-thalassemia whereby oxidative membrane remodeling of erythrocytes and iron-related hepatobiliary signaling leave distinct, measurable footprints in the serum metabolome.

A. Phosphatidylcholine biosynthesis and remodeling

Several phosphatidylcholines (PCs) and lysophosphatidylcholines (LysoPCs) were significantly elevated in β-thalassemia patients compared to healthy controls.

Phosphatidylcholines. PC(O-14:0/2:0) or 1-Tetradecyl-2-acetyl-sn-glycero-3-phosphocholine PC(0:0/18:1(9E) or 2-oleoyl-sn-glycero-3-phosphocholine, PC(17:2(9Z,12Z)/19:0) or 1,2-diacyl-sn-glycero-3-phosphocholine, PC(20:4(5Z,8Z,11Z,14Z)/18:0) and PC(18:0/20:4(8Z,11Z,14Z,17Z)) are PCs (glycerophospholipids) in which a phosphorylcholine moiety occupies a glycerol substitution site and are intermediates in the metabolic pathway of PCs. PCs and sphingomyelin are the dominant phospholipids of the erythrocyte outer membrane layer, whereas phosphatidylethanolamine and phosphatidylserine (PS) reside largely on the inner layer. This transbilayer asymmetry—maintained by flippases/floppases/scramblases—is essential for RBC deformability and for evading macrophage recognition [19]. In β-thalassemia, excess α-globin chains, labile iron, and elevated intracellular Ca²⁺ promote lipid peroxidation and scramblase activation, causing PS externalization that accelerates splenic clearance and contributes to vascular/thrombotic complications [20]. Furthermore, phosphocholine and phosphatidylcholine catabolism play a key role in ATP production and terminal erythropoiesis and the deregulation of this process may contribute to ineffective erythropoiesis in thalassemia [21].

Lysophosphatidylcholines (LysoPCs). LysoPC(18:1(11Z)) and LysoPC(18:3(9Z,12Z,15Z)) are monoacyl-glycerophospholipids bearing a phosphorylcholine headgroup. In plasma, LPCs arise primarily from PC hydrolysis by lipoprotein-associated phospholipase A₂ (Lp-PLA₂) and from lecithin:cholesterol acyltransferase (LCAT) during cholesterol esterification. Beyond their structural role, LysoPCs function as bioactive lipids: they engage lysophospholipid-responsive pathways, generally upregulating genes involved in cholesterol biosynthesis while suppressing programs of hepatic fatty-acid oxidation, and at elevated concentrations they provoke endothelial dysfunction and oxidative stress. Excess LysoPCs production can reflect increased PLA₂ flux, including Lp-PLA₂ overexpression/activation [22]. Notably, β-thalassemia patients exhibit elevated circulating Lp-PLA₂ activity [23], plausibly driven by monocyte/macrophage activation and heightened lipoprotein/erythrocyte oxidative stress. Together with our data, these observations support a mechanistic model in which erythroid and systemic redox imbalance expands the pool of oxidized PCs; PLA₂ activities (including Lp-PLA₂) cleave these substrates to generate LysoPCs that accumulate in the circulation and, in turn, amplify endothelial and inflammatory signaling.

Arachidonic-acid (AA). We detected enrichment of arachidonic-acid–linked pathways, with AA-related features among the top discriminants. AA is a PUFA esterified predominantly in phosphatidylethanolamine, phosphatidylcholine, and phosphatidylinositol pools of cellular membranes. Upon activation, PLA₂ hydrolyzes the sn-2 ester bond to liberate free AA, which is then converted by cyclooxygenases (COX-1/2), lipoxygenases (5-/12-/15-LOX), and cytochrome P450 epoxygenases into eicosanoids—potent autacoids governing inflammation, vascular tone, platelet function, and leukocyte trafficking [24].

In our cohort, the up-regulation of arachidonic-acid methyl ester and related signals indicates accelerated PUFA turnover and a shift toward pro-inflammatory oxylipin biosynthesis, consistent with elevated oxidative stress and sterile inflammation. These findings align with prior untargeted metabolomics in β-thalassemia that reported broad perturbations in lipid and energy pathways and a chronic pro-inflammatory milieu (9). In thalassemia patiens, chronic hemolysis and iron-driven reactive oxygen species (ROS) production may both increase PLA₂-dependent AA liberation and promote non-enzymatic PUFA peroxidation yielding nonclassic eicosanoids, such as isoprostanes, that further amplify vascular inflammation [23,25]. Although our untargeted approach cannot resolve individual eicosanoids, the combined signatures—elevated LysoPCs, perturbation of PUFA-rich PCs, and pathway enrichment—indicate increased flow of the PLA₂/AA axis.

B. Bile acid biosynthesis

We observed up-regulation of conjugated primary bile acids - glycocholic acid (GCA), glycolchenodeoxyglycocholic acid (GCDCA) and and of the C27 intermediate (25R)-3α,7α-dihydroxy-5β-cholestan-27-oyl taurine alongside the secondary species glycoursodeoxycholic acid (GUDCA). Primary bile acids are synthesized from cholesterol in hepatocytes and conjugated to glycine/taurine before biliary secretion; most are reabsorbed in the ileum, but a fraction undergoes microbial deconjugation and 7α-dehydroxylation/epimerization to form secondary bile acids e.g., ursodeoxycholic acid (UDCA), GUDCA, deoxycholic acid (DCA) [26]. The elevation of bilirubin and primary bile acids in our patients suggests coupled effects of hemolysis-driven heme turnover and hepatic stress. The increased levels of bile acids may have several systemic effects such as deterioration of hemolysis [28].

Functionally, excess bile acids can worsen red-cell fragility and exert system-wide metabolic and immunologic effects via the nuclear receptor FXR (farnesoid X receptor) and TGR5/GPBAR1 (Takeda G protein-coupled receptor 5/G protein-coupled bile acid receptor 1) receptor—pathways that regulate FGF19 (Fibroblast growth factor 19) signaling, lipid and glucose metabolism, energy homeostasis, and macrophage-mediated inflammation [28]. Notably, iron overload has been shown to suppress hepatocellular FXR signaling and disrupt bile-acid homeostasis, increasing circulating biosynthetic intermediates [28]—coherent with our detection of the taurine-conjugated C27 intermediate above.

The presence of GUDCA implicates active microbiome remodeling of the bile-acid pool: UDCA/GUDCA arise from bacterial epimerization of chenodeoxycholic acid followed by hepatic re-conjugation. This aligns with emerging β-thalassemia microbiome–metabolome data linking intestinal bilirubin/stercobilin turnover and specific taxa (e.g., Cupriavidus metallidurans) to anemia severity and hepatic/splenic iron burden [16]. Collectively, these results support a gut–liver axis in β-thalassemia wherein iron-related hepatobiliary signaling and microbiome-driven bile-acid transformations intersect with lipid/oxidative pathways, complementing the membrane-lipid remodeling signature seen elsewhere in our dataset.

5. Conclusions

Untargeted UHPLC–Orbitrap MS metabolomics of β-thalassemia serum samples revealed a dual metabolic signature: (i) disrupted glycerophospholipid metabolism with elevated PC/LysoPC species and (ii) altered bile acid biosynthesis marked by increased conjugated bile acids and a taurine-linked C27 intermediate. These findings suggest that oxidative stress from iron overload and hemolysis drives membrane remodeling and inflammation, while hepatobiliary and gut-derived signals reshape the bile acid pool. Notably, these metabolic alterations were consistent across TD and NTD patients. The results propose potential biomarkers (LysoPCs, AA-derivatives, bile acids) and therapeutic targets (oxidative stress, PLA₂ pathways, FXR/TGR5 signaling, microbiome modulation). Larger longitudinal studies integrating metabolomics with clinical outcomes are needed for biomarker validation and clinical translation.

Limitations of the study include the modest sample size; the wide age range; the potential confounding by transfusion timing, chelation, diet, and microbiota not fully controlled; and the inherent constraints of untargeted annotation. Finally, our cross-sectional design may under-power subgroup contrasts (TDT vs. NTDT) and cannot resolve causality.