Clinical information was obtained from the attending physician from 09 to 31, 2023. All meaningful findings were coded to the standard terminology HPO before being analyzed using the in-house phenotype-genotype pipeline IDeRare. ¹⁷

Cryopreserved DNA samples were isolated from blood samples taken from the patient and the patient’s parents using the Geneaid™ DNA Isolation Kit (Blood) according to the manufacturer’s protocol on May 20, 2021. Sequencing was conducted at the Beijing Genomics Institute (BGI) Genomics Laboratory, China, using DNA Nanoball Sequencing (DNBSeq) and Agilent V6 kit with 100x coverage on March 26, 2023, with Material Transfer Agreement (MTA) approved by the Ministry of Health Indonesia (HK.07.01/H/2443/2023). The resulting sequencing files in FASTQ format were analyzed, annotated, screened, and prioritized using the IDeRare pipeline. ¹⁷ Novel VUS was defined as rare variants with minor allele frequency < 0.05 with a lack of population statistics information or functional effects of mutations in the population. Population databases used to confirm the novelty of this variant were the gnomAD3.1 database ¹⁸ ( https://gnomad.broadinstitute.org/) and dbSNP. ¹⁹ Functional prediction of the mutations was performed using the dbNSFP 4.4a ²⁰ annotation database.

In this study, the wild-type (WT) GBE sequence was obtained from NCBI (GenPept #NP_000149.4) and edited to create novel R198T, known pathogenic L224P, and Y329S amino acid sequences. WT, L224P, Y329S, and R198T were modeled using AlphaFold2, ²¹ validated with MolProbity, ²² and superimposed using PyMol. ^{23– 25} Structural analysis of mutations was compared to the normal GBE protein residue domain, which was grouped into the N-terminal alpha helix (43–75), carbohydrate binding molecule 48 (CBM48, 76–183), central catalytic core (184–600), and C-terminal amylase-like domain (601–702) with a ligand binding site involving surface proteins from both CBM48 and the central catalytic core (62–63, 91–93, 118–121, 333–335). ^{6, 26}

Identification of potential changes in the primary structure (chemical bonding, hydrophobicity, and steric hindrances), secondary structure (protein misfolding), and tertiary structure (change in protein surface) between novel R198T mutations and normal (WT) or known-pathogenic mutations (L224P, Y329S) was performed using ExPASy, ²⁷ PSIPRED, ²⁸ and visualized using PyMol. Prediction analysis of proteasome degradation through internal and external proteasome systems was performed using the proteasome cleavage prediction server (PCPS), ²⁹ with the ubiquitination tendency of lysine (Lys, K) in the ubiquitin-proteasome system (UPS) pathway ³⁰ analyzed using RUBI. ³¹

Maltoheptose (Glc7) [Glygen #G41288IQ] ⁶ was used to model the binding site of the glucose polymer in the GBE. Potential binding sites and cavities were predicted using the ICMPocketFinder. Protein-ligand docking was performed using Autodock Vina with grid coordinate centers C _x, C _y, C _z (−13.514, 7.876, 27.993) and box dimensions of s _x, s _y, s _z (40, 40,40), centered near the ligand-binding site of normal GBE (GBE WT). The highest-rank complexes with the lowest Gibbs free energy were visualized using the BIOVIA Discovery Studio Visualizer.

HEK293T cells were grown in Dulbeccos’s modified Eagle’s medium (Gibco) with 10% fetal bovine serum (FBS), 100 U/mL penicillin, and 100 mg/mL streptomycin at 37 °C with 5% CO ₂. The cells were transiently transfected with plasmids from Synbio Technologies (USA) encoding eGFP-vector/empty vector (EV), WT, L224P, Y329S, and R198T using Lipofectamine™ 3000 Transfection Reagent (Invitrogen™ # L3000008) according to the manufacturer’s protocol. The cells were grown overnight (20 h) in basal medium before harvesting. The expression of β-actin and GBE from triplicate-independent experiments was analyzed by western blotting.

The cells were lysed in RIPA buffer supplemented with a protease inhibitor cocktail (Roche #05892970001) at a ratio of 5 × 10 ⁶ cells/1 mL of buffer. The lysis process was performed on ice for 15 min before being centrifuged at 12,000 rpm for 5 min at 4 °C. The supernatant was collected and quantified using a DC Protein Assay (Biorad #500–0116) and read on an ELISA reader at a wavelength of 655 nm. Equal mass of protein (10 μg) supernatant was boiled in Laemmli Sample Buffer for 5 min. The samples were separated using SDS-PAGE and transferred to a polyvinylidene difluoride membrane (Merck Millipore #IPVH00010) at 100 V for 1 h. The membranes were blocked in 5% milk (BIORAD #1706404) in TBS-Tween for 1 h before being incubated with 1:1000 primary antibodies rabbit anti-GBE (Thermo #PA5–26515) and 1:2500 rabbit anti- β-actin (Thermo #PA5–85271) overnight at 4 °C. The membranes were then washed three times with TBST, labeled with 1:10000 horseradish peroxidase (HRP)-conjugated secondary antibodies goat anti-rabbit (Thermo #31460), and detected by enhanced chemiluminescence (ECL) (Elabscience #E-IR-R307) using ImageQuant™ LAS 4000. The band signal was analyzed using ImageJ, and the GBE (~72 kDa) signal was normalized to the β-actin (~42 kDa) signal.

GraphPad Prism 9 was used for the analysis of in vitro results. The normality test was done using the Shapiro-Wilk test, and Analysis of Variance (ANOVA) was used to compare band intensities between mutants and wild type. Post hoc analysis was conducted using Tukey’s multiple comparisons for multiple comparisons between each group. Statistical significance was defined as P < 0.05.

Patient A was a 3-year-9-month-old Javanese boy referred from a regional referral hospital, Dr. Kariadi Hospital to the Nutrition and Metabolic Department of Child Health, Dr. Cipto Mangunkusumo National Referral Hospital (RSCM), Jakarta on March 24, 2021 because of hepatosplenomegaly since the age of 2 years with a history of weight faltering, refractory anemia, and generalized weakness. He was second-born to an asymptomatic non-consanguineous 37-year-old father and 32-year-old mother. The patient was delivered full-term by C-section with normal birth weight and length. His elder brother passed away at 2 years of age with an enlarged abdomen, which was the cause of death due to suspected thalassemia. His maternal aunts were known to have unexplained death before the age of 5 years ( Figure 1). Before his visit, he had consumed ursodeoxycholic acid (2 × 75 mg), folic acid (1 × 1 mg), and albumin (3 × 1000 mg) daily.

Physical examination at the clinic revealed anemic conjunctiva, enlarged liver 4 cm below the costal arch, spleen in the Schuffner 4 position, and enlarged abdominal circumference of 57 cm with no specific dysmorphism. Anthropometric examination showed mid upper stature (14.3 cm, Z-score − 2 to −1 SD), and normal stature (94.5 cm, Z-score − 2 to −1 SD). Laboratory workup showed anemia (10.5 g/dL), thrombocytopenia (79 x 103/μL), increased transaminase enzyme (ALT 136 U/L, AST 308 U/L), GGT of 82 U/L, cholestasis (total|direct|indirect bilirubin 3.17|1.74|1.43 mg/dL), hypoalbuminemia (3.13 g/dL), PT/APTT: 1.46x/1.5x, and normal fasting blood glucose (68 mg/dL). Follow-up blood metabolomic result: mild hyperlactatemia (3.0 mmol/L), AFP 26.09 ng/mL, metabolic alkalosis (pH 7.541, pCO2 29.6 mmHg, pO2 51 mmHg, HCO3 25.6 mmol/L, Total CO2 26.5 mmol/L, BE 4.2 mmol/L, O2 saturation 89.7%), mild hyponatremia (sodium 132 mEq/L, potassium 3.9 mEq/L, chloride 102.8 mEq/L, calcium 8.1 mg/dL, calcium ion 1.17 mmol/L, phosphate 4.3 mg/dL, magnesium 1.73 mg/dL), normal ammonia (85 μg/dL), and low HDL (Total cholesterol 122 mg/dL, HDL 26 mg/dL, LDL 85 mg/dL, Triglyceride 67 mg/dL).

Before referral, the patient was treated in the hematooncology department of a regional hospital from his first admission at the age of 2 years, before finally being referred to the nutrition and metabolic department at 3-year-2-month-old. The patient underwent numerous initial and follow-up diagnostic procedures in regional hospitals, including routine hematological, chemical, and immunological workup; hemoglobin electrophoresis; bone marrow biopsy; abdominal US and MSCT; bone survey; echocardiography; liver biopsy; hepatobiliary scan; and electromyography. Liver biopsy showed a red purplish color at hepatocyte cytoplasm inclusions on periodic acid-Schiff (PAS), with progression of liver cirrhosis. Bone marrow biopsy showed the existence of foam cells supporting lysosomal storage disease. A nerve conduction study showed a proximal lower motor neuron disease origin. Bone scans only showed osteopenia, echocardiography showed mild mitral regurgitation, abdominal US and MSCT supported hepatosplenomegaly, and hepatobiliary scans excluded biliary atresia. Several differential diagnoses were considered and ruled out upon diagnostic workup, including thalassemia, autoimmune hepatitis, common cause of liver injury or infective hepatitis, hepatoblastoma, biliary atresia, myelodysplastic syndrome, Gaucher disease, and Niemann Pick type C. Proband genetic screening was conducted after the exclusion of all potential differential diagnoses when the patient was 3-years and 5-months-old. All notable clinical findings were coded in HPO (Supplementary Table 1).

A novel GBE1mutation (NM_000158.4:c.593G > C:p. R198T) was identified by proband and trio sequencing of the heterozygous parent (ACMG PP4, Figure 2A). Variant-based prioritization conducted using IDeRare ¹⁷ detected no other known pathogenic mutations in GBE1. The frequency of this variant in the population has not been recorded in gnomAD3.1(+) or reported in NCBI dbSNP (ACMG PM2, Table 1). GBE1 R198 is a conserved region of GBE1. Annotation analysis using SIFT4G, PolyPhen2, MutationTaster, PROVEAN, MetaSVM, MetaLR, MetaRNN, and MVP consistently supported the damaging nature of the R198T mutation ( Table 1). The REVEL ensemble score of 0.873 supported moderate pathogenicity ¹⁶ in silico (ACMG PP3, Table 1). In addition, based on August 09, 2024, variant data, most missense GBE1 variants (44 out of 50) were classified as pathogenic or likely pathogenic (ACMG PP2). ³²

Molecular modeling was performed to explore the effect of the R198T mutation on the overall structure of the novel GBE1 R198T mutant compared with normal GBE1WT, pathogenic L224P, and Y329S mutants. Base substitution (G replaced by C) resulted in a coding change from arginine (Arg, R) to threonine (Thr, T) at the 198 ^th amino acid residue. The R198T mutation is located in the central catalytic core domain of GBE, which plays an important role in catalyzing glycogen branch formation, similar to L224P and Y329S ( Figure 2C).

The AlphaFold models of GBE1 R198T, L224P, Y329S, and WT were validated and superimposed. The overall modeled structure was well-modeled, with Ramachandran favoring ≥95% and MolProbity Score ≥ 66th percentile (Supplementary Table 2), with Ramachandran Plot attached in Supplementary Figure 1–3. The superimposition RMSD values of all mutants were < 2 Å, indicating minimal deviation of the overall structure of the SNV mutant compared to that of the GBE WT ( Figure 2D-F, Supplementary Table 2).

The significant change in the isoelectric point (pI) of GBE1 R198T was caused by the loss of positively charged arginine, which was replaced by threonine (neutral charge). GBE1 Y329S experienced a greater reduction in molecular weight than the other mutants due to the loss of the tyrosine phenol ring. The aliphatic index of GBE1 L224P was decreased by the loss of leucine, which has an aliphatic side chain, and its replacement with proline, which has a cyclic side chain. The Y329S mutation, compared to other mutants, exhibited hydrophobicity properties that were close to those of GBE1 WT (GRAVY Y32S -0.374 vs. WT −0.375), reflecting the preserved hydrophobicity of the GBE1 Y329S protein. All GBE mutants were categorized as stable in vitro (Supplementary Table 3).

The results of the internal proteasome and immunoproteasome cleavage site analysis showed R198T degradation change(s) at the 198 ^th residue, L224P at the 223 ^rd and 224 ^th, and Y329S at the 328 ^th and 329 ^th (Supplementary Table 4). The ubiquitination prediction results showed a slight decrease and increase in the ubiquitination tendency of R198T and L224P, respectively, compared to WT. Y329S had the same ubiquitination probability as WT, indicating a relatively similar ubiquitination pattern between Y329S and WT (Supplementary Table 5).

The GBE1 R198T mutation led to the disruption of polar contacts and loss of steric hindrance formed by G231 and E591. The polar contact with C234 is shortened (from 2.5 Å to 2.4 Å). Å N233 increased from two contact points (2.3 Å and 2.1 Å) to three contact points (2.3 Å, 2.1 Å, 2.3 Å), as shown in Figure 3A-B. The disruption of the bond with G231 and E591 caused disarrangement of the GBE catalytic core protein domain with the alpha helix, CBM, and ligand-binding site ( Figure 4A) and secondary structure misfolding (Supplementary Table 6) due to the loss of the strong bond between arginine and negatively charged amino acids such as E591. The short side chain of T198 prevented it from bonding with G231, and stronger polar contacts with C234 and N233 changed the protein surface curvature near T198 ( Figure 4G-H).

Similarly, the GBE1 Y329S mutation caused a loss of the hydrophobic phenol core, the disruption of polar contact with S290, and the formation of a new bond with F327 (2.7 Å) ( Figure 3E-F), resulting in disarranged protein domain ( Figure 4B) and replacement of hydrophobic and polar buried cavities into single protrusion polar cavity ( Figure 4E-F). Meanwhile, GBE1 L224P mutation disrupted polar contacts and increased steric hindrance with F219 due to the cyclic ring of P224 ( Figure 3E-F). Under normal conditions (WT), L224 was not an amino acid located beneath the surface ( Figure 4C-D) and the L224P mutation did not cause the disarrangement of the GBE protein domain because the existence of polar bonds with F219, I227, and K228 was retained ( Figure 3C-D).

The docking predictions were validated by comparing the rank 1 docking conformation of GBE1 WT with a previous study conducted by Froese et al. (2015). The docking results showed Glc7 contact GBE1 WT at amino acids 62–63, 91–93, 111–122, 332, and 336, similar to the findings reported by Froese et al. (2015) (Fig. S5). GBE1 L224P (Fig. S6) and R198T ( Figure 5) mutations caused changes in the Glc7 ligand-binding sites. GBE1 Y329S (Fig. S7) mutants had relatively the same ligand-binding site as GBE1 WT involving amino acids both in the CBM48 domain (62–63, 68–69, 71–72, 76, 85, 89, 91, 93–94, 116, 119–121) and central catalytic core domain (332–333, 336) despite its disarranged surface ( Figure 4B).

The predicted binding site for the GBE1 R198T-Glc7 ligand ( Figure 5) is located around R80 (red circle). The GBE1 R198T-Glc7 ligand binds through hydrogen bonds with TYR73, LYS162, ASP179, GLU181, SER183, TYR184, GLU185, PHE186, and GLU347, with weak non-covalent carbon-hydrogen interactions on HIS79, ARG80, PRO180, HIS182, ARG190, TRP344, LEU346, and GLU348, and donor-donor interactions on GLU74 and ARG350. Both L224P and R198T exhibited changes in ligand-binding sites located closer to the catalytic core, with the Glc7 ligand conforming in a circular manner on the cavity wall, in contrast to WT and Y329S, which adopt relatively ligand-receptor binding. The complete information on the available cavities on the R198T, L224P, and Y329S surfaces and protein-ligand Gibbs energies is presented in Supplementary Table 7.

Based on visual observations of the treatment group ( Figure 6A), the β-actin band appeared relatively consistent across the columns, indicating a balanced load of protein across each lane. The GBE band showed lower signal intensity for the GBE1 R198T variant than GBE1 WT, indicating that protein expression had decreased. When compared to previously known pathogenic variants, GBE1 R198T protein expression was stronger than GBE1 L224P, but weaker than GBE1 Y329S. Normalization of GBE signals resulted in band intensity of GBE1 WT: 124880 ± 53946, L224P: 10406 ± 3092, Y329S: 24903 ± 11855, and R198T: 13878 ± 7605, respectively. The Shapiro-Wilk test showed normality across three independent experimental sets ( P > 0.05), and the ANOVA test resulted in P = 0.0029. Subsequent multiple comparison tests across groups showed that GBE1 WT signal was significantly higher than all mutants: R198T ( P = 0.0054), L224P ( P = 0.0045), and Y329S ( P = 0.0100) ( Figure 6B, Supplementary Table 8), which supported the ACMG PS3 criteria for the R198T novel mutation. However, there were no significant differences between the R198T vs. L224P mutants ( P = 0.9986), R198T vs. Y329S ( P = 0.9606), or L224P vs. Y329S ( P = 0.9175) (Supplementary Table 6).

While phenotype-genotype analysis was able to match GBE1 R198T as a potential disease-causing variant, it was still classified as a VUS because this mutation only fulfilled one moderate and three supporting criteria (PM2, PP2, PP3, and PP4). The ACMG pathogenic mutation must be established by at least two moderate pathogenic and two supporting criteria (2 PM and ≥ 2 PP) or one strong pathogenic with two supporting criteria (1 PS and ≥ 2 PP criteria). Due to a lack of fulfilling moderate criteria, a single strong pathogenic criterion, such as PS3, was chosen. The protein expression study showed a lower GBE1 R198T novel variant band intensity compared to the GBE1 WT, fulfilling the PS3 criterion, which allowed us to reclassify the GBE1 R198T variant as likely pathogenic (1 PS criterion and ≥ 2 PP criteria).

Our multi-omics diagnostic experiment for GSD IV successfully demonstrated the pathogenicity and pathomechanism of GBE1 R198T. The experimental workflow was conducted step-by-step, starting with phenotype analysis, followed by genotype analysis and protein expression functional studies. As there is a lack of multiomics modalities in Indonesia and possibly in other developing countries, we proposed a modification of multiomics diagnostic workflows by prioritizing phenotype and genotype analysis before advanced omics studies. The modified workflow emphasized the importance of completing phenotype information and basic laboratory or radiological work-up within 1–2 weeks, before deciding the need for genotype analysis. Given the chance of encountering novel variants in rare disease cases, broad variant screening through WES was preferred over WGS, considering cost efficiency. Trio WES was proposed to identify patients’ VUS inheritance patterns and enhance the variant calling results. Genotype analysis can be completed within two weeks, considering that WES typically takes 5–10 days, with proband and trio analysis using the IDeRare pipeline requiring another day (~18 h).