Abstract

   Glycogen storage disease type Ia (GSDIa) is an inherited disorder of
   carbohydrate metabolism. Patients present with excessive storage of
   glycogen and fat in the liver and kidneys and are potentially at risk
   of developing long-term complications. Currently, the mainstay of
   treatment is highly tailored dietary regimens aimed at improving
   metabolic control. In the present study, to better elucidate the
   mechanisms potentially involved in the development of long-term
   complications, a mass spectrometry-based strategy was employed for an
   in-depth characterization of the serum proteomic and metabolomic
   profile of n.12 GSDIa patients. The detection of differential abundance
   of highly liver-specific circulating proteins and choline-related
   metabolites in patients provides new insights into the extent of liver
   damage and dysregulation of lipid metabolism in GSDIa. Specifically,
   the differential abundance of serum aldolase B and its positive
   correlation with traditional liver function markers supports its role
   as a potential biomarker for long-term monitoring of GSDIa liver
   injury.

   Keywords: Glycogen storage disease type Ia, Multiomics, Lipid
   metabolism, Liver injury, Serum biomarkers, ALDOB

   Subject terms: Biochemistry, Metabolomics, Proteomics

Introduction

   Glycogen storage diseases (GSDs), also referred to glycogenoses,
   constitute a group of rare, monogenic disorders caused by defects in
   glycogen synthesis or glycogen breakdown. These diseases are
   characterised by abnormal tissue deposition of glycogen^[48]1–[49]3.
   Glycogen storage disease type Ia (GSDIa; OMIM#232200) is caused by
   inactivating mutations in the G6PC1 gene encoding
   glucose-6-phosphatase-α (G6Pase-α), the catalytic subunit of the
   glucose-6-phosphatase complex. G6Pase-α plays a pivotal role in glucose
   homeostasis by catalysing the hydrolysis of intracellular
   glucose-6-phosphate (G6P) to glucose and free phosphate in the lumen of
   the endoplasmic reticulum, in the final step of gluconeogenesis and
   glycogenolysis^[50]4,[51]5. As a result, during periods of fasting,
   patients with GSDIa experience an inability to release glucose into the
   bloodstream, resulting in hypoglycaemia^[52]6. G6P is the first
   intermediate of intracellular glucose metabolism and plays a central
   role in hepatic energy metabolism, acting as a metabolic link among
   glycolysis, the pentose phosphate pathway, glycogen synthesis, de novo
   lipogenesis and the hexosamine pathway^[53]7. Consequently, GSDIa
   patients present a wide range of additional biochemical manifestations,
   such as hyperlactatemia, hyperlipidaemia, hypertriglyceridemia and
   hyperuricemia^[54]8. The tissue expression of G6Pase-α, which is
   restricted to the liver, kidneys, and, to a lesser extent, the
   intestine, determines the clinical phenotype of the disease. This
   primarily results in excessive storage of glycogen and fat in the
   tissues, leading to hepatomegaly with hepatic steatosis and kidney
   disease. Despite the recent focus on novel therapeutic approaches in
   the field^[55]9,[56]10, current treatment strategies for GSDIa remain
   largely symptomatic. The primary treatment strategy for patients with
   GSDIa involves adherence to a highly personalised dietary regimen
   comprising frequent feeding, uncooked corn starch (UCCS), and/or
   continuous nocturnal gastric drip-feeding (CNGDF). This approach has
   contributed to improvements in prognosis and long-term outcomes in
   recent decades^[57]6,[58]8,[59]11. However, over time, despite good
   metabolic control, acute hypoglycaemia and complications such as
   hepatocellular adenomas (HCA), which may progress to hepatocellular
   carcinoma (HCC)^[60]8,[61]12,[62]13, and renal failure^[63]14,[64]15
   may occur. In this regard, elevated levels of circulating lipids have
   been observed to be a contributing factor. Indeed, an association has
   been demonstrated between elevated levels of serum triglycerides (TG)
   in childhood and the risk of developing HCA in patients with
   GSDIa^[65]16, as well as an association between hyperlipidaemia and the
   progression of renal damage^[66]17. In the present study, the
   characterisation of the serum lipidomic profile previously performed in
   a cohort of patients with GSDIa^[67]18 was complemented by metabolomic
   and proteomic analysis of the same serum sample from the same patients.
   This approach provided a more comprehensive and detailed perspective on
   the metabolic changes associated with GSDIa that cannot be detected by
   conventional analytical methods. The ultimate goal of this study was to
   gain further insight into the alterations in lipid metabolism
   previously observed in GSDIa patients and to expand the investigation
   to the identification of serum metabolites and/or proteins potentially
   associated with the disease phenotype. This approach aimed to uncover
   novel biomarkers linked to GSDIa, which could contribute to a better
   understanding of the molecular mechanisms underlying the disease and
   its long-term complications. Given the complexity of the biological
   system in both health and disease, the added value of our study was a
   broader and more comprehensive knowledge in the field of GSDIa by
   integrating data from different molecular levels.

Results

Serum proteomic profile of patients with GSDIa reveals altered liver-specific
and lipid metabolism-related proteins

   In the present study, a label-free quantitative (LFQ) proteomics
   approach was employed to investigate the GSDIa serum proteome in a
   cohort of 12 patients with GSDIa, compared to 12 age- and sex-matched
   healthy controls (HC). A total of 406 proteins were identified and
   quantified, and a detailed list of the quantified proteins is provided
   in the supplementary material (Supplementary Table 1). This includes
   the UniProtKB/Swiss-Prot ID, gene name, relative quantification in each
   sample and differential analysis results for each protein. To assess
   changes in the serum proteome of patients with GSDIa and HC, entire
   data sets were processed by multivariate and univariate statistical
   analyses. Firstly, a principal component analysis (PCA) was conducted
   to estimate the variability between and within each data set. The PCA
   plot (Fig. [68]1a) illustrates the distribution of samples along the
   first two principal components (PC1 and PC2), which explain 18.06% and
   11.44% of the total variance in the data, respectively. The plot
   demonstrates effective clustering between patients and healthy controls
   (HC), with one patient (P6) borderline between the two conditions.
   Subsequent to this, a heat map was generated following hierarchical
   clustering of the most differentially expressed serum proteins for adj.
   p < 0.01 and absolute fold change |FC|> 2 between patients and HC (Fig.
   [69]1b). Based on these 24 highly differential proteins, the dendrogram
   structure exhibited clear separation of patients in the two groups,
   with the exception of patient P6. The complete heat map of differential
   proteins according to less stringent criteria (adj. p < 0.05,
   |FC|> 1.5) is provided in Supplementary Fig. [70]1a. In order to
   provide a more comprehensive representation of the serum proteins that
   exhibited the most significant differential abundance between patients
   with GSDIa and HC, a volcano plot analysis was employed, based on the
   relative abundance, expressed as log2FC, and statistical significance
   levels, expressed as − log10(p value) (see Fig. [71]1c and
   Supplementary Fig. [72]1b). Specifically, Fig. [73]1c provides the
   expression pattern of proteins differential for adj. p < 0.05 and
   |FC|> 1.5. Furthermore, the Fig. [74]1c also provides the labelling of
   24 differential proteins for adj. p < 0.01 and |FC|> 2. Among these, n.
   15 were found significantly upregulated proteins (ALDOB, APOA5, APOC3,
   APOC4, BTD, CES1, CTSD, F10, F11, FAH, ICAM1, MUC5B, PIGR, PRG4, VNN1)
   and n. 9 were significantly downregulated (ACTA1, ACTG1, APOD, CDH13,
   IGLV5-45, PI16, PZP, SHBG, VCAM1) in the serum of patients with GSDIa
   compared to HC. The levels of these differentially abundant proteins in
   GSDIa patients and HC serum samples were expressed using individual
   scatter plots (Fig. [75]1d).

Fig. 1.

   [76]Fig. 1
   [77]Open in a new tab

   Serum proteomic profile of patients with GSDIa. (a) Principal component
   analysis (PCA). The first principal component (PC1) explains 18.06% and
   the second principal component (PC2) explains 11.44% of the total
   variance in the data between patients and HC. (b) Hierarchical
   clustering with heat map visualization of differential proteins for
   adj. p < 0.01 and |FC|> 2. Heat map colour coding represents relative
   protein abundance: red and blue represent increased and decreased
   levels of each protein in patients compared to HC. (c) Volcano plot of
   differential proteins for adj. p < 0.05 and |FC|> 1.5. Labels only for
   significantly different proteins with adj. p < 0.01 and |FC|> 2. Blue
   and red dots refer to significantly decreased and increased proteins,
   respectively, in patients compared to HC. (d) Scatter plots of all
   differential proteins with adj. p < 0.01 and |FC|> 2 in terms of
   log2(LFQ) normalized values. Each dot is annotated with patient number
   and sex (♂ represents male, ♀ represents female). Missing values,
   imputed for the differential analysis are not shown. The corresponding
   tabular data are presented in Supplementary Table [78]1.

   In order to elucidate the pathophysiological mechanisms underlying
   GSDIa, the set of proteins identified as dysregulated in patients was
   subjected to further investigation. To this end, an evaluation of
   tissue gene expression and an assessment of functional enrichment
   analysis were conducted. This investigation yielded two key results.
   Firstly, bioinformatic analysis of the tissue gene expression profiles
   revealed an enrichment of proteins (ALDOB, APOC3, APOC4, APOA5, CTSD,
   CES1, ACTG1, FAH, VNN1, F11, PRG4 and F10) with predominantly hepatic
   expression (Fig. [79]2a). This finding contributed to the
   identification of the possible association between the altered levels
   of these proteins found in the patients’ serum and an organ-specific
   dysfunction. A heat map illustrating the tissue transcript expression
   profile of the differential proteins is shown in Fig. [80]2a. The
   yellow/blue colour scale for the transcript-expression-pattern (TPM,
   Transcripts Per Million) score provides information on both the
   distribution of transcript levels among tissues and the putative
   tissue-specific relative expression. To further emphasize this result,
   liver-specific transcriptomic data were extracted and visualized (Fig.
   [81]2b), highlighting the hepatic expression of the proteins included
   in the heatmap of Fig. [82]2a. In order to assess the association
   between the altered serum proteome profile and liver damage in GSDIa,
   correlation analyses were then performed between the proteins
   identified by GTEx analysis with liver TPM > 5. (ALDOB, APOC3, APOC4,
   APOA5, CTSD, CES1, ACTG1, FAH, VNN1, F11, PRG4, F10, SHBG, BTD, PIGR,
   ICAM1, PZP) and traditional biochemical liver monitoring parameters
   (ALT and AST). Both ALDOB and FAH showed a positive linear correlation
   with ALT and AST levels; however, ALDOB was selected for further
   consideration due to its higher liver expression (TPM) compared to FAH
   (Fig. [83]2b). Notably, ALDOB showed significant correlations with ALT
   (p < 0.001) (Fig. [84]2c) and AST (p < 0.01) (Fig. [85]2d) in patients
   compared to healthy controls. Furthermore, ALDOB levels showed a strong
   positive correlation with serum uric acid (p < 0.05) (Fig. [86]2e), a
   parameter commonly used to monitor kidney function and gout risk in
   GSDIa. Importantly, ALT, AST and uric acid levels were all
   significantly elevated in patients compared to controls (Fig.
   [87]2f–h), supporting the clinical relevance of these findings.

Fig. 2.

   [88]Fig. 2
   [89]Open in a new tab

   Investigation of differential serum protein content of patients with
   GSDIa. (a) Heatmap of tissue gene expression profile of differential
   proteins in GSDIa. Heatmap color coding represents different transcript
   distribution values in different tissues, expressed as TPM. Clustering
   of liver proteins was observed. (b) Extrapolation of liver-specific
   tissue gene expression profile from (a). (c–e) Simple linear regression
   analysis between log2(LFQ) normalized ALDOB values and serum levels of
   ALT (c) and AST (d), and uric acid I in patients with GSDIa. Solid
   lines represent the result of the linear fit. Solid line colour code
   represents patient (red) and HC (grey) trend. R: Pearson’s correlation.
   P[unadj]: uncorrected p value. Missing values, imputed for the
   differential analysis are not shown. (f–h) Plots of differential serum
   ALT (f), AST (g), and uric acid (h) levels, expressed as (means ± SEM)
   in GSDIa patients versus HC. The significant difference (*p < 0.05,
   **p < 0.01, ***p < 0.001) was evaluated by performing Mann–Whitney
   comparison test. (i) Lollipop charts of the top 10 biological processes
   enriched from the input data list of differentially expressed proteins
   in the serum of patients with GSDIa. The −log10(FDR) values for each
   pathway were expressed on a color scale between orange and purple. (j)
   Plot of differential serum total cholesterol levels, expressed as
   (means ± SEM) in GSDIa patients versus HC. The significant difference
   (*p < 0.05, **p < 0.01, ***p < 0.001) was evaluated by performing
   Mann–Whitney comparison test.

   Secondly, a functional enrichment analysis of the differential protein
   list revealed a significant enrichment of processes primarily related
   to lipid metabolism using an FDR cutoff of 0.05 and a minimum pathway
   size of 5. In order to identify specific altered biological processes
   underlying the GSDIa condition, the Gene Ontology (GO) database was
   used to analyse the data. The top 10 enriched pathways are shown in
   Fig. [90]2i. Consistently with these molecular findings, Fig. [91]2j
   reports a plot showing significantly altered levels of total
   cholesterol from routine clinical biochemistry, highlighting a marked
   increase in GSDIa patients compared to healthy controls (p < 0.0001).

Serum metabolomic profile of patients with GSDIa reveals features consistent
with liver dysfunction and dysregulation of lipid metabolism

   In the present study, the serum metabolomic profile of patients with
   GSDIa was obtained by using a targeted liquid chromatography–tandem
   mass spectrometry (LC–MS/MS)-based approach, which allowed to robustly
   detect and quantify 73 metabolites. A comprehensive list of the
   identified and quantified metabolites, including metabolite names,
   abbreviations, group classifications and their raw concentrations, is
   provided in Supplementary Table 2. Changes in the serum metabolome of
   patients with GSDIa and HC were investigated using multivariate and
   univariate statistical analyses.

   To estimate the variability between and within each group of samples,
   data were analyzed by PCA. With PC1 and PC2 variances of 25.85% and
   17.35% respectively, the PCA score plot showed a clear separation
   between patients with GSDIa and HC along the first PC (Fig. [92]3a).
   Metabolites selected as significantly differentially concentrated in
   sera of patients with respect to HCs were visualized in a heat map
   after hierarchical clustering (Fig. [93]3b). With the exception of one
   healthy control (C6), the dendrogram structure showed clear clustering
   of the two patient groups. The results of the differential analysis are
   reported in Supplementary Table 2, and significantly abundant
   metabolites (adj. p < 0.01, |FC|> 1.5) were also shown using a volcano
   plot to highlight significance as a function of relative abundance of
   each metabolite, expressed as − log10 (p value) and log2 fold change,
   respectively (Fig. [94]3c). The analysis identified 14 out of 73
   metabolites as differentially abundant between the serum samples of the
   two study groups. More specifically, patients with GSDIa showed higher
   levels of beta-Ala, carnosine, choline, Glu, Ala, Asp, Pro, lactate
   (Lac), sarcosine, serotonin, xanthine, and lower levels of
   p-cresol-SO4, Phenylalanine betaine and Trimethylamine N-oxide (TMAO),
   compared to HC. The log2-transformed concentrations of the
   significantly differentially abundant metabolites are presented as
   scatter plots in Fig. [95]3d. Among these, the altered levels of
   choline, sarcosine, and TMAO—metabolites intricately linked to
   phospholipid synthesis and lipoprotein metabolism—point to a potential
   dysregulation of lipid metabolic pathways in GSDIa patients. Moreover,
   correlation analyses were conducted in order to explore potential
   associations between amino acids abundance in patients and their
   circulating levels of ALT and AST. Specifically, significant positive
   linear correlations were found between serum levels of Ala and ALT
   (p < 0.05), Ala and AST (p < 0.01), Asp and ALT (p < 0.01), Asp and AST
   (p < 0.001), Glu and ALT (p < 0.05) and Glu and AST (p < 0.05) (Fig.
   [96]4a–f), These findings suggest a strong link between elevated levels
   of these amino acids and traditional liver function parameters.

Fig. 3.

   [97]Fig. 3
   [98]Open in a new tab

   Serum metabolomics profile of patients with GSDIa. (a) Principal
   component analysis (PCA). The first principal component (PC1) explains
   25.85% and the second principal component (PC2) explains 17.35% of the
   total variance in the data between patients and HC. (b) Hierarchical
   clustering with heat map visualization of differential metabolites for
   adj. p < 0.01 and |FC|> 1.5. Heat map colour coding represents relative
   metabolite abundance: red and blue represent increased and decreased
   levels of each metabolite in patients compared to HC. (c) Volcano plot
   of differential metabolites for adj. p < 0.01 and |FC|> 1.5. Blue and
   red dots refer to significantly decreased and increased metabolites,
   respectively, in patients compared to HC. (d) Scatter plots of the
   significantly differential metabolites in terms of log2 raw
   concentration. Missing values, imputed for the differential analysis
   are not shown. The corresponding tabular data are presented in
   Supplementary Table [99]2.

Fig. 4.

   [100]Fig. 4
   [101]Open in a new tab

   Relevant clinical correlations between differential serum metabolites
   and traditional parameters of liver function in GSDIa patients. (a–f)
   Simple linear regression plots of significant correlations found
   between differential serum Ala, Asp, Glu levels and serum ALT, AST
   levels in patients with GSDIa. Solid lines represent the result of the
   linear fit. Solid line colour code represents patient (red) and HC
   (grey) trend. R: Pearson’s correlation. P[unadj]: uncorrected p value.
   Missing values, imputed for the differential analysis are not shown.

Multiomics data integration provides new insights into relationships gained
and lost in disease

   To highlight pairs of features that could be directly or indirectly
   functionally related or take part in common biological processes,
   Pearson’s correlation coefficients were calculated for each
   protein–protein, metabolite-metabolite, and protein-metabolite pair.
   Correlations were calculated separately for GSDIa patients and HCs and
   considered significant with an adjusted p value < 0.1. All significant
   correlations in either HCs or patients are shown in Supplementary Fig.
   [102]2a, b, respectively, and in Supplementary Table 3. Furthermore, to
   systematically investigate the role of GSDIa disease on these
   correlations, those gained or lost under pathologic conditions were
   reported in Supplementary Fig. [103]2c and Supplementary Table [104]3,
   and the most significant gained or lost correlations were sub-selected
   and shown in Fig. [105]5.

Fig. 5.

   [106]Fig. 5
   [107]Open in a new tab

   Multiomics correlation changes in GSDIa. (a) Selected lost or gained
   significant Pearson’s correlations in patients with respect to HCs
   (adjusted p-value < 0.1) with |corr.coeff patient – corr.coeff.
   HC|> 0.6 and at least one of the connected features significantly
   different between patients and HCs. Pink circular sectors correspond to
   proteins, green to metabolites. The width of the circular sectors is
   proportional to |corr. Coeff. Patient – corr. Coeff. HC|. The orange
   and purple edge colours indicate gained and lost correlations in
   patients. Molecules highlighted in red are significantly increased in
   patients, those highlighted in blue are decreased. (b) Scatter plot of
   z-scores of the indicated features with significant correlation
   difference between control and patient. Solid lines represent the
   result of the linear fit. Solid line colour code represents patient
   (red) and HC (grey) trend. R: Pearson’s correlation. P[unadj]:
   uncorrected p value. Missing values, imputed for the differential
   analysis are not shown.

   In HC serum samples, 180 significant correlations were found
   (Supplementary Fig. [108]2a), and, out of these, only 12 were multiomic
   (between a protein and a metabolite). In particular, negative
   correlations were found between SHBG and Asparagine, Leucine,
   Isoleucine, and Lysine, between Desmocollin 1 (DSC1) and Symmetric
   dimethylarginine, between Superoxide Dismutase 3 (SOD3) and Tryptophan
   betaine, between Suprabasin (SBSN) and Serotonin, and between
   Immunoglobulin Lambda Variable 4-69 (IGLV4-69) and
   Dehydroepiandrosterone sulfate. Positive correlations were found
   between CBP2 (a.k.a. SERPINH1) and α-, β-, and γ-Aminobutyric acid, and
   between Complement Factor H Related 3 (CFHR3) and Tryptophan betaine.
   All these multiomic correlations in HC samples were lost in GSDIa
   samples (Supplementary Fig. [109]2c). For GSDIa patient serum samples,
   140 significant correlations were found (Supplementary Fig. [110]2b),
   with only 3 of these being multiomic: between
   Alpha-1-Microglobulin/Bikunin Precursor (AMBP) and Aspartate, between
   endoplasmic reticulum chaperone BiP (HSPA5) and Homoarginine, and
   between complement factor H-related protein 2 (CFHR2) and Lactic acid,
   which were all negative correlations and not significant in HC samples
   (“gained correlations” in Supplementary Fig. [111]2c).

Discussion

   The life expectancy of patients with GSDIa has increased significantly
   over the past decades as a result of advances in dietary
   management^[112]11. Nevertheless, despite the evidence indicating that
   good metabolic control plays a role in preventing complications in
   GSDIa^[113]19, acute and chronic complications still occur in patients
   despite adherence to dietary therapy^[114]8. In this study, a broad
   metabolomic approach was applied with the aim of identifying additional
   metabolic disruptions and gaining deeper insights into the systemic
   biochemical imbalances associated with GSDIa. Serum samples from 12
   patients with GSDIa and sex- and age-matched healthy controls were
   analyzed. The metabolomic and proteomic serum profiles of GSDIa
   patients were characterised, allowing a detailed assessment of
   disease-related molecular alterations.

   The characterization of the serum proteomic profile in patients with
   GSDIa enabled the identification and quantification of 406 proteins.
   Among these, 24 were found to be significantly dysregulated compared to
   healthy controls. Enrichment analysis of these differentially expressed
   proteins revealed two main features. First, a strong involvement in
   lipid metabolism emerged, consistent with previous findings from our
   earlier study^[115]18 which characterised the serum lipidomic profile
   of the same cohort of GSDIa patients analysed here. Second, a
   predominant liver-specific expression pattern was observed, further
   supporting the central role of hepatic dysfunction in the molecular
   landscape of the disease.

   Moreover, the results of the present study are consistent with the
   presence of liver damage in GSDIa. In this regard, a more detailed
   examination of the literature revealed that the dysregulated proteins
   identified in the serum of patients with GSDIa overlap with the
   alteration profile observed in the context of liver disease. This
   finding underscores the potential for shared mechanisms underlying
   GSDIa and liver disease. In particular, it has been reported that the
   levels of eight serum proteins (ALDOB, PIGR, SHBG, APOD, PRG4, APOB,
   PZP and IGLV5-45) vary in individuals with liver diseases other than
   GSDIa.

   Among these, ALDOB, the hepatic isozyme of fructose-bisphosphate
   aldolase (EC 4.1.2.13), and PIGR, an immunoglobulin receptor involved
   in the transcytosis of IgA and IgM, have recently been proposed as
   potential circulating biomarkers of metabolic dysfunction-associated
   steatotic liver disease (MASLD) in humans^[116]20. Specifically, both
   ALDOB and PIGR have been found to be dysregulated in liver diseases,
   including MASLD, metabolic dysfunction-associated steatohepatitis(MASH)
   and alcohol-related liver disease (ALD). In detail, ALDOB was found to
   be up-regulated in plasma from both patients and a mouse model of
   MASLD^[117]21, in plasma from individuals with MASH^[118]22, in serum
   from patients with acute and chronic hepatitis^[119]23, and
   simultaneously up- and downregulated in plasma and liver from
   individuals with ALD^[120]24. PIGR was found to be upregulated in the
   plasma of both patients with MASLD ^[121]21,[122]25 and a mouse model
   of MASLD^[123]21, in the liver and plasma of patients with ALD^[124]24,
   and in the plasma of patients with MASH^[125]22 and MASH rat
   model^[126]26. Moreover, ALDOB was found to be upregulated in both the
   plasma and liver of ob/ob mice, a leptin-deficient model of obesity,
   whereas PIGR was only upregulated in their plasma^[127]27. It is
   noteworthy that, in accordance with the findings of this study, ALDOB
   was also identified as being upregulated in the livers of two distinct
   GSDIa mouse models^[128]28,[129]29. In the present study, in alignment
   with the findings of ALDOB in the context of liver disease, we found a
   linear positive correlation between circulating ALDOB levels and the
   transaminases AST and ALT, which are traditional biochemical parameters
   of liver function; furthermore, consistent with ALDOB implication in de
   novo lipogenesis in MASLD^[130]30, we also found a positive correlation
   between ALDOB and the fructose-related metabolite uric acid (Fig.
   [131]2e), suggesting increased fructose metabolism in GSDIa to provide
   precursor for gluconeogenesis and de novo lipogenesis. Our previous
   work demonstrated a specific altered serum lipid signature in the same
   cohort of patients^[132]18.

   Similarly to ALDOB and PIGR, there is evidence to suggest that also the
   sex hormone binding globulin (SHBG), which was found to be
   downregulated in GSDIa patients in the present study and in Simons et
   al.^[133]31, may also play a role in the pathophysiology of MASLD. In
   more detail, retrospective and prospective studies of very large
   cohorts of patients with MASLD^[134]32,[135]33, have associated serum
   levels of SHBG with the development and regression of MASLD, suggesting
   that SHBG is a new promising biomarker of MASLD. Indeed, SHBG is
   recognized as a hepatokine and a reliable marker of MASLD, as well as
   other conditions, including obesity, diabetes mellitus, and insulin
   resistance^[136]34. SHBG is a homodimeric glycoprotein, produced by
   hepatocytes and secreted into the blood, where it binds with high
   affinity to circulating sex steroid hormones, and regulates their
   bioavailability^[137]35. In this regard, the dysregulation of SHBG and
   the demonstration of an inverse correlation between SHBG and the
   saturated fatty acids fraction and the intrahepatic lipid content in
   individuals with GSDIa^[138]31, provide evidence supporting the
   hypothesis that there is a link between metabolic imbalance in GSDIa
   and the circulating sex hormone levels in GSDIa which has recently been
   the subject of increased attention^[139]36. Indeed, additional
   features, including delayed puberty, hypogonadotropic hypogonadism and
   polycystic ovaries are documented in GSDIa^[140]36. In regard to this
   topic, the findings of a comprehensive systematic review^[141]37 and a
   meta-analysis indicate a correlation between reduced SHBG levels and
   the presence of polycystic ovary syndrome^[142]38. It is noteworthy
   that, although dietary influence on growth and sexual development is
   evident, the imbalance in circulating levels of insulin and cortisol
   observed in hepatic GSDs, whether secondary to dietary treatment or to
   the enzyme defect, may play a pivotal role in this
   sense^[143]36,[144]39. Indeed, chronic recurrent elevations of cortisol
   in response to hypoglycemia may lead to the suppression of
   gonadotropin-releasing hormone (GnRH), luteinizing hormone (LH), and
   follicle-stimulating hormone (FSH) release, with implications for SHBG
   biosynthesis^[145]36.

   It should be noted that the present study found differentially abundant
   proteins in the serum of patients with GSDIa, including polymeric
   immunoglobulin receptor (PIGR), proteoglycan 4 (PRG4), pregnancy zone
   protein (PZP) and immunoglobulin lambda variable 5-45 (IGLV5-45). These
   proteins may also be relevant in the context of liver disease. These
   have been linked to HCC studies, and are joined by the previously
   discussed ALDOB protein. In particular, ALDOB has been identified as a
   potential prognostic biomarker for overall survival in HCC in a recent
   study using big data analysis^[146]40. PIGR, which was found to be
   upregulated in the serum of patients with GSDIa in the current study,
   is also being investigated for its ability to exert oncogenic functions
   through the activation of the ribosomal pathway in HCC^[147]41 and
   metastasis^[148]42. PRG4 is a secreted proteoglycan protein
   predominantly expressed by chondrocytes, which lubricate and protect
   the surfaces of articular cartilage^[149]43. Its mRNA expression in
   liver cancer tissue was recently shown to correlate with longer
   survival in HCC patients likely due to its ability to impair HCC cell
   migration and act as a tumor suppressor in tumors expressing the major
   PGR4 receptor, CD44^[150]44. In the case of PZP, an estrogen-related
   protein that normally increases during pregnancy, there is evidence
   that its downregulation in HCC tissue may be associated with immune
   infiltrates and poor prognosis^[151]45. Furthermore, IGLV5-45, the
   fifth region of the variable domain of immunoglobulin lambda light
   chains, which was found to be downregulated in our patients, was found
   to be upregulated in study of immunodeficiency herpes virus/hepatitis B
   co-infected patients with HCC^[152]46.

   This study identified differentially abundant proteins in the serum of
   patients with GSDIa, suggesting possible liver damage, as supported by
   the aforementioned considerations. Further analysis revealed that some
   of these differentially abundant proteins have a role in lipid
   metabolism, which is consistent with our previous lipidomic studies in
   the same cohort of patients^[153]18. Proteins that are differentially
   abundant in patients with GSDIa and implicated in lipid metabolism
   include APOD, APOA5, APOC3, APOC4, BTD, CES1, VNN1, and the above
   discussed SHBG, PIGR, PRG4, and ALDOB.

   Specifically, APOA5, APOC3, APOC4, BTD, VNN1, ALDOB and PIGR were all
   found to be upregulated in the plasma proteome of ob/ob mice, whereas
   ALDOB, APOA5 and VNN1 were also found to be upregulated in both plasma
   and liver^[154]27. In particular, the alterations observed in the
   apolipoproteins (APOD, APOA5, APOC3, and APOC4), in conjunction with
   those in the cholesteryl ester hydrolase (CES1), which is predominantly
   expressed in the liver^[155]47, provide robust evidence for lipid
   metabolism disturbances in GSDIa, as corroborated by the lipidomics
   results^[156]18. With regard to APOD, which was identified as being
   downregulated in the present study in GSDIa, it has been demonstrated
   that it is predominantly present in high-density lipoproteins (HDL) and
   capable of promoting VLDL catabolism^[157]48. In particular, increased
   APOD production has been demonstrated to significantly reduce plasma TG
   levels in mice^[158]48, while its deficiency has been shown to result
   in elevated serum TG and insulin levels in mice^[159]49. The
   downregulation of APOD in our patients is therefore consistent with
   these findings and with the results of an altered serum lipidome
   profile in patients with GSDIa^[160]18 and delayed very-low-density
   lipoprotein catabolism in GSDIa^[161]50,[162]51. Finally, also
   biotinidase (BTD) and PRG4 seem to be implicated in lipid and glucose
   metabolism perturbations. In detail and consistent with BTD result
   obtained, BTD has been proposed as a potential marker of perturbed
   glucose and lipid metabolism in GSDs^[163]2,[164]52,[165]53. In line
   with the overabundance of PRG4 found in the serum of GSDIa patients in
   this study and the metabolic disturbances typical of GSDIa discussed,
   it has recently been suggested that PRG4 contributes to weight gain,
   dyslipidaemia, and insulin resistance^[166]43.

   The characterization of the serum metabolomic profile of patients with
   GSDIa allowed the identification and quantification of 73 metabolites
   out of a total of 80 present in the targeted metabolomics panel. The
   application of univariate and multivariate statistical analysis
   revealed differential abundance of 14 metabolites in the serum of
   patients. In particular, patients with GSDIa exhibited elevated levels
   of beta-Ala, carnosine, choline, Glu, Ala, Asp, Pro, Lac, sarcosine,
   serotonin, xanthine, and diminished levels of p-cresol-SO,
   Phenylalanine betaine and TMAO in comparison to HC. The elevated serum
   concentrations of Lac are indicative of the hyperlactatemia observed in
   patients with GSDIa. This condition is believed to arise from an
   augmented intracellular glycolytic rate, driven by the accumulation of
   G6P, as previously suggested^[167]54. With regard to the levels of Glu,
   Ala, Asp, and Pro found in patients, these amino acids can act as
   gluconeogenic amino acids in a variety of processes involving
   transamination to pyruvate (Ala) and through shunting into the TCA
   cycle by transamination to α-ketoglutarate (Glu), to oxaloacetate
   (Asp), and by catabolism to α-ketoglutarate (Pro). In particular,
   transamination reactions are catalysed by ALT and AST, whose serum
   levels appear to be elevated in patients included in this study.
   Furthermore, linear positive correlations were identified between serum
   ALT and Ala, Glu, and Asp, as well as between serum AST and Asp and
   Ala. These findings indicate a robust correlation between these amino
   acids and aminotransferase levels. In particular, the comparison with
   HC revealed a greater and more significant difference in ALT levels
   than in AST levels. As ALT assays measure the combined activity of both
   ALT1 and ALT2 isoenzymes, and since ALT1 is predominantly
   liver-specific^[168]55, the elevated ALT levels observed likely reflect
   liver damage in patients with GSDIa. It is noteworthy that serum
   concentrations of AST and ALT are elevated in GSDIa at the time of
   diagnosis and often return to normal or near-normal levels with
   appropriate treatment^[169]8. Nevertheless, and despite good compliance
   with dietary treatment, our results show elevated circulating AST and
   ALT levels in patients with GSDIa, in line with the findings of Farah
   et al^[170]56 and Rossi et al^[171]54 in patients with GSDIa. These
   results provide further evidence that dietary intervention alone is not
   sufficient to prevent liver damage.

   To provide a more detailed explanation, elevated serum levels of Ala
   and Asp, when considered as products of transamination reactions, may
   indicate an increased glycolytic rate (associated with elevated
   pyruvate levels) and TCA cycle rate (associated with elevated
   oxaloacetate levels) or an increased need for transport to the liver to
   supply the gluconeogenic pathway. In particular, there are
   findings^[172]54,[173]56 indicating the presence of an increased
   glycolytic rate and TCA cycle overload in patients with GSDIa. In sum,
   it can be speculated that in patients with GSDIa, despite optimal
   glycemic control, there is a positive regulation of gluconeogenesis by
   both TCA cycle overload (mediated by Asp, Glu, and Pro) and
   transamination to pyruvate (mediated by Ala) in response to subclinical
   hypoglycemia. It should also be noted that, in line with data
   suggesting TCA overload in patients^[174]54,[175]57, higher levels of
   fumarylacetoacetase (FAH) were detected by serum proteomics in patients
   with GSDIa in the present study.

   The elevated serum concentrations of xanthine observed in this study
   are consistent with the untargeted plasma metabolomics analysis
   conducted by Mathis et al.^[176]56 in 14 patients with GSDI, including
   11 with GSDIa. This corroborates the hypothesis of altered purine
   metabolism in GSDIa^[177]56.

   Furthermore, the present metabolomic study revealed an elevation in
   serum choline and sarcosine levels, alongside a pronounced reduction in
   TMAO, in patients with GSDIa in comparison to HC. Choline is an
   essential biological molecule for all cells, required for the synthesis
   of phosphatidylcholine (PC), a structural component of
   VLDL^[178]58,[179]59. PCs participate in the SM cycle and play a key
   role in the ceramide (Cer) /sphingomyelin (SM) balance, through the
   hydrolysis of SM by sphingomyelinase and the transfer of phosphocholine
   from PC to Cer by sphingomyelin synthase^[180]60. In addition, Choline
   is also a source of labile methyl groups in the metabolism of
   homocysteine, through the conversion to betaine^[181]61. In accordance
   with this, the present study identified elevated serum concentrations
   of sarcosine and TMAO, two choline-related metabolites, and carnosine
   (and its precursor, beta-Ala). Carnosine is involved in the reaction
   catalysed by carnosine N-methyltransferase to produce anserine and
   S-adenosyl-homocysteine. Sarcosine is one of the choline-related
   metabolites^[182]62. Specifically, sarcosine is derived from choline
   through a series of chemical reactions, including the oxidation of
   choline to form betaine. Subsequently, betaine undergoes conversion to
   dimethylglycine, which in turn is converted to sarcosine. TMAO is a
   choline-related metabolite, resulting from the oxidation of
   trimethylamine by flavin monooxygenase 3 in the liver. A reduction in
   TMAO may be indicative of alterations in choline metabolism and
   dysbiosis of the microbiota in patients with GSDIa, as demonstrated
   recently^[183]63. The findings of the present study have been used to
   construct a working model of a pathomechanism affecting choline
   metabolism in GSDIa (Fig. [184]6). This model integrates the current
   metabolomic results with the results of our previous study on
   lipidomics^[185]18. The present study hypothesises that an increase in
   hepatic choline metabolism in patients with GSDIa supports the increase
   in PC levels in patients via the Kennedy pathway. This increase, in
   turn, is hypothesised to support an increase in the SM rate and Lands
   cycle. Indeed, in our previous study^[186]18 we demonstrated an
   increase in PC, DG, SM, CER, and LysoPC serum levels in the same cohort
   of patients with GSDIa.

Fig. 6.

   [187]Fig. 6
   [188]Open in a new tab

   Proposed pathomechanisms involving choline metabolism in GSDIa. The
   increase in choline and the decrease in the choline-related metabolite
   TMAO, in parallel with an increase in the choline-related metabolite
   sarcosine, and the higher levels of PC, LysoPC, DG, SM and Cer obtained
   in our previous serum lipidomics study on the same cohort of patients
   with GSDIa, suggest an increased rate of the Kennedy pathway and,
   simultaneously, of the SM and Lands cycles in GSDIa. In this context,
   the role of PC in the Cer/SM balance is crucial: phosphocholine is
   transferred from PC to Cer via SMS and subsequently released as
   diacylglycerols. This process is coupled to the hydrolysis of SM to Cer
   by sphingomyelinase. In addition, PC promotes the increase of LysoPC in
   the Lands cycle through its hydrolysis mediated by PLA2. Choline is
   also a methyl group donor in homocysteine metabolism through its
   conversion to betaine and its metabolism to dimethylglycine. The
   increase in serum sarcosine levels from dimethylglycine can be
   explained as compensation for the increased levels of choline converted
   to PC in the Kennedy pathway. Metabolites detected and identified as
   elevated are shown in red; those detected and identified as decreased
   are shown in blue. The presence of an asterisk refers to metabolites
   that were found to be differentially abundant in our previous
   lipidomics study18 performed on the same cohort of GSDIa patients.
   Metabolites that were measured and found not to be differentially
   abundant are indicated by a yellow colouring. Abbreviations. ALDH7A1,
   betaine aldehyde dehydrogenase; BHMT, betaine-homocysteine
   S-methyltransferase; CCT, phosphocholine cytidylyltransferase; CHDH,
   choline dehydrogenase; CK, choline kinase; CPT,
   cholinephosphotransferase; DMGDH, dimethylglycine dehydrogenase; FMOs,
   flavin monooxygenase enzymes; LPCAT, Lysophosphatidylcholine
   acyltransferase; PLA2, Phospholipase A2; SMS, sphingomyelin synthase.

   The following study assessed, through a multi-omics approach, potential
   correlations between metabolites and proteins in order to identify
   meaningful associations related to the disease. This strategy allowed
   the identification of molecular relationships that may arise from
   shared pathophysiological mechanisms, even without a direct link. It is
   important to note that correlation does not imply causation. For
   instance, two features may appear correlated because they are both
   indirectly influenced by a third factor, rather than being directly
   related to each other. There is no direct evidence in the literature
   linking many of the identified correlated features and follow-up
   research is required to unravel the underlying mechanisms.
   Nevertheless, correlation analysis—both at the single- and multi-omics
   levels—can help identify disease-associated processes, even when the
   molecules involved are differently affected, as observed in the case of
   sarcosine and hypoxanthine. The present study revealed a significant
   increase in serum sarcosine concentration in patients, whereas
   hypoxanthine levels did not differ between groups. However, patient
   samples exhibited a higher concentration of xanthine (a metabolite of
   hypoxanthine) than HC subjects. As previously discussed, sarcosine is a
   secondary metabolite of choline, an important component of hepatic
   lipid metabolism, and xanthine oxidase, the enzyme responsible for the
   oxidation of hypoxanthine to xanthine and subsequently to uric acid, is
   predominantly expressed in the liver. Consequently, these observations
   imply that the observed correlation between sarcosine and hypoxanthine
   may be attributable to the concomitant presence of altered (albeit
   distinct) processes associated with liver disease. It is noteworthy
   that this positive relationship between sarcosine and hypoxanthine is
   only observed in the context of the disease; in HCs, no discernible
   relationship between these two molecules is observed. This observation
   can be attributed to the disease-specific nature of this relationship.

Limitations of the study

   Despite the strengths outlined above, this study has limitations mainly
   due to the rarity of GSDIa, which resulted in a small sample size. This
   limitation consequently reduced the statistical power and hindered the
   ability to conduct subgroup analyses, including sex-specific
   evaluations and assessments related to pharmacological treatments and
   dietary interventions. While the limited cohort size did not allow for
   a dedicated sex-specific analysis, we performed reanalysis including
   sex as a covariate for key markers such as PZP and SHBG, confirming
   that the results remained significant. Nonetheless, we acknowledge the
   lack of comprehensive sex-specific assessments as a limitation that
   should be addressed in future studies with larger cohorts. Although
   detailed stratification based on treatment regimens was not feasible,
   individual results were interpreted, when possible, considering the
   clinical management of each patient. For example, while allopurinol—a
   xanthine oxidase inhibitor commonly used to prevent gout in GSDIa
   patients—may have confounding effects, elevated serum xanthine levels
   were observed even in four of six patients not receiving this treatment
   (P1, P2, P3, P11), suggesting that other factors contribute to these
   metabolic alterations. Nutritional status represents an additional
   potential confounder influencing metabolic profiles, warranting
   cautious interpretation of the results. In particular, the findings
   regarding SHBG levels in our patients may have been influenced by
   overnutrition^[189]64,[190]65, which is often a consequence of the
   specific dietary regimens implemented to manage hypoglycemia in GSDIa.
   Moreover, previous studies have shown that plasma choline levels show
   minimal variation in healthy individuals^[191]66, and fasting-mimicking
   diets have been reported to reduce plasma TMAO levels^[192]67
   supporting the view that the metabolic changes observed in our cohort
   are primarily disease-related. Nevertheless, these aspects underscore
   the importance of thoroughly evaluating dietary influences in future
   research to better define the role of the identified biomarkers.
   Furthermore, the lack of longitudinal data limits our ability to assess
   the stability and intra-individual variability of these biomarkers over
   time, which is crucial for confirming their applicability in ongoing
   disease monitoring.

Conclusions

   The objective of this study was to provide a foundation for the
   advancement of current research priorities in the domain of hepatic
   glycogenosis. The investigation was designed to add to the existing
   body of knowledge concerning metabolic alterations in GSDIa, which are
   not discernible through conventional monitoring analyses. Complementing
   the findings of our previous lipidomic analysis^[193]18 the
   comprehensive assessment provided by this multiomics approach provided
   new insights into the pivotal role of the liver in the pathophysiology
   of GSDIa and a detailed assessment of lipid metabolism dysfunction in
   GSDIa. In this regard, the serum metabolomics findings on choline
   metabolism support our previous lipidomics data and thus contribute to
   the conceptualization of a novel working model of choline metabolism
   impairment in GSDIa. In addition, the positive linear correlation of
   ALDOB with patient ALT and AST supports further research to investigate
   its role as a novel potential biomarker of liver injury. Finally, the
   positive relationship between ALDOB and patient serum uric acid levels
   supports the potential role of ALDOB in providing substrates for de
   novo lipogenesis in GSDIa. In conclusion, this study underscores the
   translational potential of multi-omics approaches to advance the
   monitoring of GSDIa. The biomarkers identified here offer promising
   avenues for improving the management and early detection of long-term
   complications associated with the disease. Further validation in
   longitudinal cohorts will be essential to confirm their clinical
   utility and integration into future guidelines.

Methods

Human participants

   This study was conducted in accordance with the principles of the
   Declaration of Helsinki. All studies were performed after informed
   consent which was obtained from adult participants or the parents of
   infants.

   For patients followed up at the University Medical Center Groningen
   (Groningen, The Netherlands), the local Medical Ethical Committee
   indicated that the Medical Research Involving Human Subjects Act was
   not applicable and official study approval by the Medical Ethical
   Committee (METc) was not required according to protocol n. METc
   2019/119. The experimental protocols were approved by the METc at the
   University Medical Center Groningen (Groningen, The Netherlands), in
   accordance with the relevant guidelines and regulations included in the
   n. METc 2019/119 protocol. For patients followed up at the University
   of Naples Federico II (Naples, Italy), the study protocol was approved
   by the Ethics Committee of the University of Naples Federico II
   (Naples, Italy), n. 77/21. The experimental protocols were approved by
   the Ethics Committee at the University of Naples Federico II (Naples,
   Italy), in accordance with the relevant guidelines and regulations
   included in n. 77/21 protocol.

   Participants were recruited over a 12 month period and the diagnosis of
   GSDIa was based on mutation analysis of the G6PC1 gene. The study
   included 12 patients with GSDIa (8 males and 4 females, median age
   16.5 ± 9 years, age range 5–34 years) and 12 age- and sex-matched HC (8
   males and 4 females, median age 20.3 ± 9.4 years, age range
   5–31.5 years). Among the patients, one was aged < 6 years, five were
   aged 10–16 years, and six were aged > 17 years. Among the HC, one was
   aged < 6 years, five were aged 10–16 years, and six were
   aged > 17 years. Table [194]1 summarises the general characteristics of
   the two study cohorts. There was no significant difference between
   patients and controls in terms of age and sex as calculated by T-test
   for age and chi-squared test for sex. Details of the clinical
   characteristics of the patients, including age, sex, genotype, dietary
   regimen, body mass index, medication and dietary supplements are shown
   in Table [195]2. The HC received no dietary or pharmacological
   treatment.

Table 1.

   General characteristics of the study cohorts.
   Parameters Subjects p value
   Patients with GSDIA Healthy controls Patients with GSDIa versus Healthy
   controls
   Size sample 12 12 –
   Age (years) 18 (5–34) 20 (5–31.5) 0.65
   Males/females 8/4 8/4 1
   Ethnicity Caucasian Caucasian –
   [196]Open in a new tab

Table 2.

   Clinical features of patients with GSDIa. ^a BMI, body mass index.
   Patients with GSDIa Age (years) Gender G6Pase variation BMI^a (kg/m^2)
   Dietary regimen Drugs—Supplements Clinical features
   P1 5 M

   p.Arg83Cys

   p.Arg83Cys
   21.0 Frequent feedings; UCCS (3.8 g/kg/day); CNGDF Vitamin D
   Intellectual disability
   P2 10 M

   p.Arg83Cys

   p.Arg83Cys
   19.7 Frequent feedings UCCS (5.2 g/kg/day) Vitamin D ^–
   P3 17 M

   p.Trp63Arg

   p.Arg83Cys
   21.4 Frequent feedings Glycosade (3.1 g/kg/day) – Previous growth
   hormone deficiency (rhGH withdrawn 6 months before study enrolment)
   P4 25.5 F

   p.Trp63Arg

   p.Arg83Cys
   24.4 Frequent feedings Glycosade (3.2 g/kg/day)

   Allopurinol

   Ramipril

   Fenofibrate
   Proteinuria
   P5 24 F

   p.Arg83Cys

   p.Arg83Cys
   22.5 Frequent feedings; UCCS (3.9 g/kg/day) ω3 fatty acids –
   P6 34 M

   p.Arg83Cys

   p.Arg83Cys
   26.7 Frequent feedings; Glycosade (4.5 g/kg/day)

   Allopurinol

   Carbamazepine
   P7 11 F

   p.Arg83Cys

   p.Arg83Cys
   22.2 Frequent feedings; UCCS (4.8 g/kg/day)

   Allopurinol

   Ramipril

   Vitamin D

   ω3 fatty acids
   Glomerular hyperfiltration
   P8 28.5 M

   p.Arg83Cys

   p.Arg83Cys
   29.6 Frequent feedings; UCCS (1.5 g/kg/day); CNGDF Oxybutynin Urinary
   incontinence
   P9 24 F

   p.Arg83Cys

   p.Arg83Cys
   26.3 Frequent feedings; UCCS (4.3 g/kg/day)

   Allopurinol

   Ramipril

   Vitamin D

   ω3 fatty acids
   Glomerular hyperfiltration
   P10 12 M

   p.Arg83Cys

   p.Arg83Cys
   19.7 Frequent feedings; CNGDF

   Allopurinol

   Vitamin D

   ω3 fatty acids
   –
   P11 10 M

   p.Gly270Val

   p.Gly270Val
   21.7 Frequent feedings; UCCS (6.3 g/kg/day) Vitamin D –
   P12 16 M

   p.Gly270Val

   p.Gly270Val
   20.2 Frequent feedings; Glycosade (2.1 g/kg/day) CNGDF

   Allopurinol

   Gemfibrozil
   Hepatocarcinoma
   [197]Open in a new tab

Serum collection

   Blood samples were collected in the pre-prandial state according to
   each patient’s fasting tolerance (3–6 h). To onooxyg bias due to short
   fasting times, HC were asked to schedule their blood sampling after the
   same fasting time as their age- and sex-matched GSDIa counterparts.
   Blood was centrifuged (10 min, 2000 × g at 4 °C) and serum was
   aliquoted into separate polypropylene tubes which were immediately
   stored at − 80 °C.

LFQ-proteomics

   S-TrapTM micro spin column (Protifi, Hutington, USA) digestion was
   performed on 50 µg of serum blood according to manufacturer’s
   instructions. Briefly, samples were reduced with 20 mM
   Tris(2-carboxyethyl)phosphine (TCEP) and alkylated with 50 mM
   chloracetamide (CAA) for 15 min at room temperature. Aqueous phosphoric
   acid was then added to a final concentration of 2.5% following by the
   addition of S-Trap binding buffer (90% aqueous methanol, 100 mM
   triethylammonium bicarbonate (TEAB), pH7.1). Mixtures were then loaded
   on S-Trap columns. Three extra washing steps were performed for
   thorough SDS elimination. Samples were digested with 2.5 µg of trypsin
   (Promega) at 47 °C for 1.5 h. After elution, peptides were vacuum
   dried.

   The tryptic peptides were resuspended in 50 µL of 2% acetonitrile (I),
   0.1% Formic acid in high-performance liquid chromatography (HPLC)-grade
   water and an amount of 200 ng was injected on a nanoElute (Bruker
   Daltonics, Germany) HPLC system coupled to a timsTOF Pro (Bruker
   Daltonics, Germany) mass spectrometer. HPLC separation (Solvent A: 0.1%
   formic acid in water; Solvent B: 0.1% formic acid in I) was carried out
   at 250 nL/min using a packed emitter column (C18, 25 cm × 75 μm 1.6 μm)
   (Ion Optics, Australia) using a gradient elution (2 to 11% solvent B
   during 19 min; 11 to 16% during 7 min; 16% to 25% during 4 min; 25% to
   80% for 3 min and finally 80% for 7 min to wash the column).
   Mass-spectrometric data were acquired using the parallel accumulation
   serial fragmentation (PASEF) acquisition method. The measurements were
   carried out over the m/z range from 100 to 1700 Th. The range of ion
   mobilities values from 0.65 to 1.4 V.s/cm^2 (1/k0). The total cycle
   time was set to 1.16 s and the number of PASEF MS/MS scans was set to
   10.

   The obtained data were analyzed using MaxQuant version 2.0.1.0 and
   searched with Andromeda search engine against the UniProtKB/Swiss-Prot
   Homo sapiens database (release 02-2021, 20,404 entries). To search
   parent mass and fragment ions, we set a mass deviation of 3 ppm and
   20 ppm respectively. The minimum peptide length was set to 7 amino
   acids and strict specificity for trypsin cleavage was required,
   allowing up to two missed cleavage sites. Carbamidomethylation (Cys)
   was set as fixed modification, whereas oxidation (Met) and N-term
   acetylation were set as variable modifications. The false discovery
   rates (FDRs) at the protein and peptide level were set to 1%. Scores
   were calculated in MaxQuant as described previously. Proteins were
   quantified according to the MaxQuant label-free algorithm using LFQ
   intensities; protein quantification was obtained using at least 1
   peptide per protein. Match between runs was allowed.

Targeted metabolomics

   The serum targeted metabolomic analysis was performed using liquid
   chromatography–tandem mass spectrometry (LC–MS/MS) based methods as
   reported by Caterino et al^[198]68,[199]69. Metabolomic analysis was
   performed according to the protocols of the MxP^® Quant 500 kit
   (Biocrates Life Sciences AG, Innsbruck, Austria) on a Triple Quad™
   5500 + System-QTRAP^® Ready (AB Sciex, Framingham, MA, USA) coupled to
   a 1260 Infinity II HPLC (Agilent, Santa Clara, CA, USA). The LC–MS/MS
   based platform was set up to quantify the concentration of the
   following molecules: amino acids (n = 20), amino acid related (n = 30),
   biogenic amines (n = 9), carboxylic acids (n = 7), hormones and related
   (n = 4), indoles and derivatives (n = 4), and alkaloids, amine oxides,
   cresols, nucleobases and related, vitamins and cofactors (n = 6).
   Specifically, 10 µL of each serum sample was pipetted onto a 96-well
   extraction plate containing wells for blanks, PBS, calibrators and
   quality controls (QC) and dried under a nitrogen stream for 30 min.
   Quantification was carried out using internal standards (ISTD) and a
   7-point calibration curve (Cal 1 to Cal 7). Three human plasma samples
   spiked with different concentrations of reference analytes (QC1-3) were
   analyzed as quality control, according to manufacturer’s protocol. Each
   sample was then incubated with 50 µL of 5% phenyl isothiocyanate (PITC)
   in ethanol/water/pyridine (ratio 1/1/1, v/v/v) for 60 min to derivatise
   amino acids and biogenic amines and then dried again under a nitrogen
   stream for 60 min. Metabolites were extracted with 300 µL of 5 mM
   ammonium acetate in methanol in a shaker (for 30 min, 450 rpm) and
   eluted into a new 96-well plate by centrifugation.

   An aliquot of 150 µL of each extract was diluted with an equal volume
   of HPLC grade water and the LC–MS/MS measurements were performed to
   target and quantify metabolites by multiple reaction monitoring
   (MRM)-based method. Liquid chromatography was performed on a C18 column
   (Zorbax eclipse XDB 3 × 100 mm, 3.5 µm) using acetonitrile and water
   with 0.2% formic acid as running solvent. For the LC analysis, the
   injection volume was 5 µL at an initial flow rate of 0.05 mL/min using
   the follow gradients: 100% A, 0–0,25 min; 88% A, 0,25–1 min; 82,5% A,
   1–3 min; 50% A, 3–4,5 min; 0% A, 4,5–7,5 min for LC1 part coupled with
   positive ion mode acquisition and 100% A, 0–0,25 min; 75% A,
   0,25–0,50 min; 50% A, 0,5–3 min; 25% A, 3–4 min; 0% A, 4,5–7,5 min for
   LC2 part coupled with negative ion mode acquisition. The autosampler
   was cooled at 10 °C. The ion source was operated in positive ion mode
   with the following parameters: spray voltage 55 kV, temperature 450 °C,
   GS1 20, GS2 40, CUR 30, CAD 8. The ion source was operated in negative
   ion mode with the following parameters: spray voltage 55 kV,
   temperature 450 °C, GS1 20, GS2 40, CUR 30, CAD 8.

   Data were generated using Analyst software v.1.7.1 (AB Sciex,
   Framingham, MA, USA) and further processed to calculate metabolite
   concentrations using MetIDQ™ Oxygen software (Biocrates Life Sciences
   AG, Innsbruck, Austria).

Proteomics bioinformatic analysis

   Bioinformatic analysis and results visualization was performed in R v.
   4.3. Reverse proteins, contaminants, and proteins only identified by
   site were filtered out. Proteins that were not quantified in at least 8
   samples of at least one condition (HC or patient) were also filtered
   out. The remaining 406 proteins were further analyzed. Data were log2
   scaled. Missing values were imputed by impSeqRob^[200]70 with alpha
   parameter equal to 0.95 (less than 5% outlying proteins). Sample-wise
   normalization was performed by median to obtain the final normalized
   log2 (LFQ) data used for differential analysis. Z-scores for heat map
   visualization were obtained by autoscaling. PCA was performed by a
   singular value decomposition of the centered data matrix using function
   prcomp. Differential analysis between HC and patient conditions was
   performed by SAM method^[201]71 within the R package samr v. 3.0, using
   t-test and 1000 permutations. Significantly differential proteins were
   identified with a mixed criterion based on adjusted p value (adj. p)
   and fold change (FC), as indicated in specific figure captions.
   Significant features that had missing (imputed) data were also
   reanalyzed using sex as a covariate to confirm significance, excluding
   an unbalance of sex samples between controls and patients. Full data
   are presented in Supplementary Table 1. Hierarchical clustering with
   heat map data visualization was performed by package pheatmap v.
   1.0.12, using Euclidean distance and complete clustering method.
   Enrichment analyses were performed using GTEx portal and ShinyGO v.
   0.81. Tissue-specific gene expression data were retrieved from the GTEx
   (Genotype-Tissue Expression) database, which provides transcriptomic
   profiles (TPM values) across a wide range of human tissues. These data
   were used to evaluate the tissue distribution of differentially
   expressed proteins and to generate heat maps of expression patterns.
   Pathway enrichment analysis was carried out using ShinyGO v. 0.81 to
   identify significantly overrepresented biological processes.

Metabolomics bioinformatic analysis

   Bioinformatic analysis and results visualization was performed in R v.
   4.3. Metabolite nomenclature reported in Supplementary information was
   provided in the targeted metabolomic kit (Biocrates Life Sciences AG).
   Metabolites that were not quantified in at least 8 samples of at least
   one condition (HC or patient) were filtered out. The remaining 73
   metabolites were further analyzed. Data were log2 scaled. Missing
   values were imputed by impSeqRob^[202]70 with alpha parameter equal to
   0.95 (less than 5% outlying proteins). Z-scores for heat map
   visualization were obtained by autoscaling. The PCA was performed by a
   singular value decomposition of the centered data matrix using function
   prcomp. Differential analysis between control and patient conditions
   was performed by SAM method^[203]71 within the R package samr v. 3.0,
   using t-test and 1000 permutations. Significantly differential
   metabolites were identified with a mixed criterion based on adj. p and
   FC, as indicated in specific figure captions. Full data are presented
   in Supplementary Table 2. Hierarchical clustering with heat map data
   visualization was performed by function pheatmap v. 1.0.12, using
   Euclidean distance and complete clustering method.

Multiomics bioinformatic analysis

   Multi-omic data analysis was performed integrating proteomic and
   metabolomic data in terms of z-scores. Pearson’s correlation
   coefficients between features in HC and patient samples were calculated
   separately, using R package Correlation v. 0.8.4, with a
   Benjamini–Hochberg p value correction for controlling false discovery
   rate. Result visualization was performed by R package Circlize v.
   0.4.16, selecting only significant correlations according to an
   adjusted p value < 0.1. Significant correlations are reported in
   Supplementary Table 3. Pairs of features whose correlation changed
   significance between patients and HCs were identified and also reported
   in Supplementary Table 3. Lost correlations in patients are
   correlations which are significant in HC (adjusted p value < 0.1) and
   non-significant in patients (adjusted p value > 0.1), vice versa holds
   for gained correlations. The results in Fig. [204]5 only report the
   most significant relationships, sub-selected considering an absolute
   difference of correlation coefficients between patient and HC samples
   larger than 0.6, and constrained to have at least one significantly
   differential feature in the lost or gained correlation. Scatter plots
   were used to display the z-scores of all non-imputed values of pairs of
   features with significant changes in correlation between HC and
   patient. Linear models were fitted to HC and patient data separately.

Data availability statement

   Data supporting the findings of this study are available in the paper
   and Supplemental information.

Supplementary Information

   [205]Supplementary Information 1.^ (2.1MB, pdf)
   [206]Supplementary Information 2.^ (189.6KB, xlsx)

Acknowledgements