Integrated multi-omics reveals anaplerotic rewiring in methylmalonyl-CoA mutase deficiency

Abstract

Methylmalonic aciduria (MMA) is an inborn error of metabolism with multiple monogenic causes and a poorly understood pathogenesis, leading to the absence of effective causal treatments. Here we employ multi-layered omics profiling combined with biochemical and clinical features of individuals with MMA to reveal a molecular diagnosis for 177 out of 210 (84%) cases, the majority (148) of whom display pathogenic variants in methylmalonyl-CoA mutase (MMUT). Stratification of these data layers by disease severity shows dysregulation of the tricarboxylic acid cycle and its replenishment (anaplerosis) by glutamine. The relevance of these disturbances is evidenced by multi-organ metabolomics of a hemizygous Mmut mouse model as well as through identification of physical interactions between MMUT and glutamine anaplerotic enzymes. Using stable-isotope tracing, we find that treatment with dimethyl-oxoglutarate restores deficient tricarboxylic acid cycling. Our work highlights glutamine anaplerosis as a potential therapeutic intervention point in MMA.

Fig. 1. Multi-faceted omics view enabled a molecular diagnosis in 84% of individuals.

a, Study overview with a depiction of the propionate pathway, including its precursors and the pathways catalyzed by MMUT. b, MMUT enzyme activity per study sub-cohort (MMUT-deficient, n = 150; other MMA, n = 60; unaffected, n = 3); P values were calculated by Wilcoxon rank test, two-sided. c, Lollipop plot of all pathogenic variants found on the MMUT gene. d, Proportions of variant types as identified on the MMUT gene. e,f, Transcript and protein levels of MMUT by study sub-cohorts (number of samples for the transcript and protein plot, respectively: MMUT-deficient, n = 143/150; other MMA, n = 59/60; unaffected, n = 19/20). g, Gene ranks according to P values as calculated by gene-wise Welch’s t-test (two-sided) in the proteomics and transcriptomics data. h, Lollipop plot of pathogenic variants identified in ACSF3. i, Proportions of affected genes identified in the whole cohort. P values were calculated by the Wilcoxon rank test. Boxplot elements represent center line, median; box limits, upper and lower quartiles; whiskers, 1.5× interquartile range; and points, outliers in e and f. Dots in b,e,f represent biologically independent samples. Source data

Fig. 2. Phenomics analysis reveals two main surrogate markers of disease severity (CSS and PI+ activity).

a, Correlation matrix of all continuous numeric and discrete phenotype variables. b, Number of phenotypic traits according to five phenotype subcategories. c, Panel of selected phenotypic traits and their overall strength of representing the entirety of the phenomics dataset (here termed disease severity) as assessed by linear modeling after log transformation. Each point represents the result of linear regression against one other phenotypic variable with the effect size (ES) on the y axis and the resulting Benjamini–Hochberg adjusted P value (two-sided) on the x axis. The horizontal curved line indicates the density of data points as distributed along the x axis. The vertical dashed line indicates the threshold of significance (P < 0.05). d, Linear regression results of the PI+ activity variable compared to the rest of the phenotypic variables; P values calculated as in c. Source data

Fig. 3. Untargeted integration of omics data layers highlights the TCA cycle and associated pathways as well as oxidative phosphorylation gene sets to be dysregulated in MMA.

a, Gene set enrichment test using the multi-omics factor analysis tool (MOFA). Benjamini–Hochberg adjusted P values (two-sided) are shown. b, Detailed feature statistics of the top enriched gene sets following MOFA in the proteomics data; P values calculated as in a. c, Linear discriminant model (split to assign training and test data, 0.5) of transcripts separates MMUT-deficient from control driven by MMUT and other genes related to the TCA cycle. d, Gene set enrichment analysis based on ES ranking derived from differential expression analysis (also Fig. 4b); P values were calculated with the fgsea R package. e, Breeding scheme of Mmut-deficient mice. f, Untargeted metabolomics in mouse tissues and body fluids, depicting boxplots for methylmalonic acid and metabolite set enrichment analysis based on the complete metabolomics dataset. Significantly changing ions between mutant and control conditions were identified using a two-sample t-test. Pathway analysis was performed using an annotated ion list ranked by P value significance. Pathway enrichments were calculated using KEGG metabolic pathway definitions and a hypergeometric test. g, RNA-seq on mouse brain tissue. Boxplots of the relative Mmut transcript abundance and gene set enrichment analysis following DESeq2 analysis; dot size represents the number of genes per set. P values in f and g (boxplot) are calculated by Wilcoxon rank test, two-sided. Boxplot elements represent center line, median; box limits, upper and lower quartiles; and whiskers, 1.5× interquartile range. Source data

Fig. 4. Transcript–protein and protein–protein correlation analyses reveal coordinated relationships between MMUT and TCA genes and proteins but not their isoforms.

a, Circos plot depicting raw fold changes (FC) of transcripts and proteins, effects sizes (ES) derived from differential expression analysis, transcript–protein correlations (rho) and correlative relationships of the MMUT protein to TCA proteins and their corresponding isoforms. b, Q-Q and volcano plots illustrate the results of the differential expression analysis based on a linear mixed modeling approach applied to the proteomics and transcriptomics data, restricted to enzymes (or their encoding genes) localized in the mitochondria; calculated P values were two-sided. c, Histograms of Spearman correlations across 4,318 transcript–protein pairs (left) and 221 samples (right). d, Scatter plot of Spearman correlations in MMUT-deficient against control. Euclidean distance from the diagonal is calculated based on the formula (MMUT-deficient correlation−control correlation)/sqrt(2); calculated P values were two-sided.

Fig. 5. Polar metabolomics and glutamine tracing studies in CRISPR/Cas9 KO 293T cells and primary patient fibroblasts highlight differential glutamine anaplerosis.

a, Volcano plot depicting differentially expressed metabolites. Highlighted are those particularly relevant to this study. b, Schematic depiction of the TCA cycle and relevant anaplerotic reactions. The color code indicates dysregulations at the metabolite and protein levels; gray metabolites were not detected. c, Pool sizes of metabolites in control and CRISPR/Cas9 KO 293T cells (error bars represent s.d., centered around the mean). d, Schematic representation of labeling of TCA cycle and associated metabolites derived from labeled glutamine via anaplerosis. Circles represent carbon atoms. e, Relative abundance of isotopologs of TCA cycle metabolites after glutamine labeling. f, Ratios of M + 5 and M + 4 citrate isotopologs. g, Total pool sizes of TCA cycle metabolites under different treatment conditions. Oxoglut., dimethyl-oxoglutarate. h, Levels of propionylcarnitine in primary patient fibroblasts under treatment. i, Ratios of M + 5 and M + 4 citrate isotopologs under treatment in primary fibroblasts; P values were calculated by two-sided Wilcoxon rank test. Boxplot elements represent center line, median; box limits, upper and lower quartiles; and whiskers, 1.5× interquartile range. For experiments in 293T cells (c,e–g), n = 3 biologically independent samples (WT), n = 2 (MMUT-KO), n = 2 (DLST-KO) over two independent experiments were measured. For experiments in patient fibroblasts (h,i), n = 4 biologically independent samples per group were measured. Source data

Fig. 6. MMUT interacts physically with GLUD1, DLST and GOT2 as demonstrated by FLAG-tag pull-down.

a, Outline of experimental and control groups indicating which protein was used with a FLAG-tag in a cross-linking affinity purification experiment coupled to subsequent analysis of the pull-down samples by mass spectrometry. b, Venn diagram of proteins pulled down by the different FLAG-tagged proteins. c, ANOVA P values (two-sided) of all proteins pulled down by MMUT. Gray dots indicate proteins with mitochondrial location, according to UniProt. d, Interaction network of significantly enriched proteins (ANOVA P value <0.05, two-sided). Thicker connector indicates lower P value. Blue, proteins with FLAG-tags used in pull-down; red, significantly pulled-down proteins as indicated in red in b that do not share any peptides with the negative controls; no color, significantly enriched proteins, but not exclusively pulled down by FLAG-tagged proteins. e, Western blot of IP of FLAG-tagged MMUT probing for DLST. Data are representative of three independent experiments. Source data

Extended Data Fig. 1. Historic context of sample collection and quality control measurements of multi-omics data.

a, Histogram of fibroblast samples binned into their year at collection and waffle chart illustrating the different sample groups indicated by the color code. b, Violin plots illustrating average number of read per sample (n = 229). HQ, high quality. Line plot (one line per sample) indicating genome coverage as quantitatively summarized in the below table. c, Boxplots indicating Phred quality scores at different number of cycles. Each sample (n = 221) underwent 75 cycles that were binned (see x axis) to display the Phred scores per bin (bin 1−2, 1206 scores; 3−5, 1809; 6−10, 3015; 11−20, 6030; 21−50, 18090; > 50, 15075). d, Proteomics quality control illustrated by the number of detected proteins in relation to the percentage of overlapping proteins (top panel). Using a threshold of > 50% overlapping proteins, the variation coefficients of n = 2850 proteins are displayed for each sample group. Boxplot elements represent center line, median; box limits, upper and lower quartiles; whiskers, 1.5x interquartile range; points, outliers. Violin plots depict the distribution of the data using vertical density curves.

Extended Data Fig. 2. Biochemical assessment of MMUT activity and propionate incorporation activity supports diagnosis of affected individuals.

a, Scatter plot of maximal, that is supplemented with adenosylcobalamin (AdoCbl) or hydroxocobalamin (OHCbl), activity of the MMUT enzyme and the propionate incorporation assay. b, Boxplots of MMUT enzyme activity with and without AdoCbl supplementation measured in biologically independent fibroblast samples (mut⁰, n = 119; mut^-, n = 29; other MMA, n = 46; unaffected, n = 3). c, Boxplots of propionate incorporation activity with and without OHCbl supplementation measured in biologically independent fibroblast samples (mut⁰, n = 120; mut^-, n = 30; other MMA, n = 60; unaffected, n = 9). d, Copy number variants illustrated by read counts of specific locations of the MMUT gene for three specific samples. e, Scatter plots of MMUT transcript and protein levels of the MMUT-deficient samples. Samples are indicated by dots and are grouped according to the underlying bi-allelic genetic variation type of the MMUT gene (Number of samples for the transcript and protein plot, respectively: deletion/deletion, n = 1/1; missense/missense, n = 63/65; missense/splicing, n = 6/6; missense/truncating, n = 26/30; splicing/splicing, n = 6/6; truncating/splicing, n = 3/3; truncating/truncating, n = 36/37). f, Regression plots of MMUT transcript and protein levels versus MMUT enzyme and propionate incorporation activity. g, Same as f but excluding samples with truncating/truncating, splicing/splicing and truncating/splicing MMUT allele combinations. h, Fold change of all transcripts and proteins, respectively, when comparing the MMUT-deficient group versus the rest of the samples. Gene names are ranked according to the negative base 10 logarithm of the fold change. All linear regressions are calculated according to the Pearson method, P values two-sided; bands indicate 95% confidence level interval. Boxplot elements represent center line, median; box limits, upper and lower quartiles; whiskers, 1.5x interquartile range; points, outliers.

Extended Data Fig. 3. Expression outlier analysis reveals causative genes in specific disease samples.

Expression rank plots for a, ACSF3, b, SUCLA2, c, MMAA and d, MMAB and Z-score volcano plots for specific samples, applying the OUTRIDER R package.

Extended Data Fig. 4. The clinical severity score and propionate incorporation activity are associated with several phenotypic traits.

a, Proportional bar plots of the presence of absence of clinical parameters in relation to the clinical severity score. b, Age at onset in relation to the clinical severity score. c, Linear regression of various relationships of propionate incorporation activity to continuous phenotypic variables (Pearson method, two-sided, pairwise comparisons). d, Comparison of propionate incorporation activity with discrete phenotype variables (P values by t-test, two-sided, pairwise comparisons; number of samples per boxplot are indicated above the boxplot). Boxplot elements represent center line, median; box limits, upper and lower quartiles; whiskers, 1.5x interquartile range. Each dot represents a sample/patient.

Extended Data Fig. 5. Transcript–protein and protein–protein correlation analysis illustrates coordinated regulation of MMUT with most TCA transcripts and proteins.

a, Histograms of Pearson correlations across 4318 transcript–protein pairs (top) and 221 samples (bottom). b, Scatter plot of the strongest positive and negative transcript–protein correlations (ranked by average of Pearson correlation coefficient in the MMUT-deficient and control group). c, Transcript–protein Pearson correlation plots of selected TCA cycle related genes. d, Spearman correlation of the MMUT protein versus a selection of TCA cycle and related proteins and their isoforms illustrated in a chord plot for control (left), MMUT-deficient samples (middle) and the difference of the two former plots (right); thickness of the links indicates nominal value of the correlation coefficient. e, Scatter plot of MMUT protein versus both isoforms of aconitase (ACO1 and ACO2) with linear regression by Pearson correlation. f, Chord plot illustrating all correlative relationships of TCA cycle and related transcripts or proteins (g); thickness of the links indicates nominal value of the Spearman correlation coefficient. All P values calculated two-sided.

Extended Data Fig. 6. Metabolomics investigation of a subset of patient cell lines.

a, Density plot illustrating a model for OGDH and GLUD1 with all fibroblast samples ranked according to three sample groups. b, Boxplots of total ion current assessed in the two experimental groups ‘control’ and ‘MMUT-deficient’ primary fibroblasts (n = 6 in each group); technical replicates are collapsed to represent one dot per cell line. c, TCA metabolites as measured by untargeted polar metabolomics; n = 6 biological replicates. d, Levels of metabolites involved in the two enzymatic steps catalyzed by two oxoacid dehydrogenase complexes (OGDC and OADC) and their proximal reactions; OGDH protein in green indicates its downregulation as detected in the proteotyping dataset; n = 6 biological replicates. Boxplot elements represent center line, median; box limits, upper and lower quartiles; whiskers, 1.5x interquartile range. All P values are calculated by Wilcoxon rank test, two-sided.

Extended Data Fig. 7. Treatment of primary patient/control fibroblasts and 293T cells.

a, Relative incorporation into different isotopologues of succinate without treatment of 293T cells. b, Relative abundance of isotopologues of TCA cycle metabolites after glutamine labeling. c, Total pool sizes of TCA cycle metabolites under different treatment conditions (Oxoglut., dimethyl-oxoglutarate) in primary fibroblasts; n = 4 biologically independent samples per group were measured. d, Total ion count for citrate and malate (e) is displayed for different isotopologues under citrate and malate treatment, respectively. For each boxplot representing results from 293 T cells (c, e, f, g), n = 3 biologically independent samples (WT), n = 2 (MMUT-KO), n = 2 (DLST-KO) over 2 independent experiments were measured. For experiments in patient fibroblasts (h, i), n = 4 biologically independent samples per group were measured; P values calculated by Wilcoxon rank test, two-sided. Boxplot elements represent center line, median; box limits, upper and lower quartiles; whiskers, 1.5x interquartile range; dots in a, d, and e are outliers, in c individual samples.

Extended Data Fig. 8. Fractional labeling pattern derived from glutamine in 293T and primary fibroblast cells upon treatment.

Cells (a, 293T. b, primary patient and control fibroblasts) treated with the compounds indicated above the plots (DM-Oxog., dimethyl-oxoglutarate).

Extended Data Fig. 9. MMUT-flag is enzymatically active and pulls down other propionate pathway proteins using immunoprecipitation.

a, MMUT enzyme activity using a radio-labeled substrate in 293T cells (backgrounds (WT or MMUT-KO) and transfected constructs are indicated in the color key; data points indicate means of n = 3 independent experiments; error bars indicate SD, centered around the mean; P values calculated by t-test, two-sided. b, Western blots of over-expressed flag-tagged MMUT, MMAB, MMAA, and MCEE, but not VLCAD, ACO2, and empty vector (EV). c, Detection of endogenous MMUT and MMAB following cross-linking immunoprecipitation. Panels shown are representative of at least n = 3 independent experiments. Source data

Extended Data Fig. 10. Quantitative pull-down results following affinity purification mass spectrometry and confirmation of DLST pull-down by MMUT-flag by immunoprecipitation.

a, Numbers indicate count of proteins per intersection following affinity purification mass spectrometry; Venn regions are labeled with the names of flag-tag bait proteins; UpSet plot illustrates intersection and set sizes. b, Immunoprecipitation of flag-tagged MMUT and MCEE probing for DLST in 293T WT and MMUT-KO cell lines. Panels shown are representative of at least n = 3 independent experiments. Source data