Clinical Variability and Genotype–Phenotype Correlation in Spanish Patients with Type 1 Gaucher Disease: A Focus on Non-c.[1226A>G]; [1448T>C] Genotypes

TL;DR

This study explores how different genetic variations in Spanish Gaucher disease patients affect their symptoms, showing that some genotypes lead to more severe health issues.

Contribution

The study identifies specific GBA1 genotypes beyond the common N370S; L444P combination that correlate with more severe clinical outcomes in Gaucher disease.

Findings

Patients with N370S combined with non-L444P variants had higher rates of advanced bone disease.

Non-N370S; L444P genotypes were strongly linked to increased risk of Parkinsonism.

Genotype stratification improves prediction of skeletal and neurological complications in GD1.

Abstract

The clinical heterogeneity of type 1 Gaucher disease (GD1) underscores the limited correlation between the GBA1 genotype and phenotype. This study examined GD1 patients from the Spanish Gaucher Disease Registry carrying heterozygous GBA1 genotypes distinct from NM_000157: c.[1226A>G](N370S); [1448T>C](L444P). Among 374 patients with GD1, 195 (52.1%) had alternative heterozygous combinations, including variants corresponding to severe (37.9%) or moderate (42.1%) mutation, whereas only 20% patients harbored mild variants—all of them in combination with N370S. Descriptive statistics and predictive models based on logistic regression and decision trees were applied. Patients carrying N370S with a different L444P variant showed significantly higher rates of advanced bone disease (59.9%) compared to those with homozygous N370S (38.3%) or N370S; L444P (41.0%) (p = 0.002). Decision tree analysis identified the bone marrow burden score (S-MRI) as the strongest predictor of osteopenia/osteoporosis at diagnosis. Genotype also emerged as a key discriminator for Parkinson’s disease: patients with non-N370S; L444P genotypes showed a markedly higher likelihood of developing Parkinsonism. Overall, GD1 patients with genotypes other than N370S; L444P present more severe phenotypes, particularly with greater skeletal involvement and neurological complications. These findings highlight the importance of genotype stratification and predictive modeling in improving risk assessment and clinical management in GD1.

Figures

Full text

1. Introduction

Type 1 Gaucher disease (GD1, OMIM #230800) is a rare lysosomal storage disorder characterized by a wide spectrum of clinical manifestations, ranging from asymptomatic presentations to severe multi-organ involvement. It is caused by biallelic pathogenic variants in the GBA1 gene (MIM*606463), which encodes the lysosomal enzyme β-glucocerebrosidase [1]. Deficiency of this enzyme leads to the accumulation of glucocerebroside in macrophages, resulting in diverse systemic complications [2].

Despite substantial research, the correlation between GBA1 genotypes and clinical phenotypes remains incomplete, limiting our ability to predict disease severity and progression based solely on genetic information [3]. Numerous studies have reported intrafamilial variability among individuals with the same genotype, including siblings, underscoring the unpredictable nature of clinical outcomes [4]. Similarly, large registries and case series have documented significant phenotypic diversity among patients with identical GBA1 variants [5].

In general, the c.1226A>G variant is associated with milder disease, particularly in homozygous individuals, whereas homozygosity for the c.1448T>C variant is linked to neuronopathic forms of Gaucher disease [6]. In Spain, the most common GBA1 genotype is the heterozygous association c.1226A>G; however, other, less common genotypes are also present and exhibit considerable phenotypic variability [7]. The relationship between the severity of enzyme dysfunction and clinical expression continues to be a subject of interest in predicting individual disease trajectories. However, the extensive allelic heterogeneity and the presence of rare variants hinder the establishment of robust genotype–phenotype correlations [8].

This complexity suggests that additional genetic, epigenetic, and environmental modifiers likely contribute to the observed clinical heterogeneity in GD1. Several studies have explored the influence of other genetic loci [9], epigenetic mechanisms, and biomarkers identified through proteomic, lipidomic, and metabolomic analyses to better understand phenotypic expression [10,11,12]. Furthermore, the genotyping of GBA1 is complicated by the presence of a highly homologous pseudogene, which increases the risk of partial recombination and can result in misclassification of alleles, particularly when using next-generation sequencing (NGS) platforms [13].

Recent advances in machine learning and complex network analysis offer promising tools for deciphering the intricate relationships between genetic, clinical, and biochemical variables in GD1. These approaches enable the analysis of large, multidimensional datasets and may improve diagnostic precision, guide treatment strategies, and enhance prognostic predictions. Their application in rare diseases such as GD1 is especially valuable as they can provide novel insights into the mechanisms driving disease variability [14,15,16,17].

In this context, the present study aims to assess the clinical characteristics at diagnosis and during follow-up in a cohort of Spanish GD1 patients with heterozygous GBA1 genotypes other than the common NM_000157: c.[1226A>G]; [1448T>C]. By focusing on this subgroup, we aim to improve our understanding of the genotype–phenotype correlations and contribute to the broader knowledge of factors influencing clinical variability in GD1.

2. Results

2.1. Patients Genotype and Characteristics

A total of 358 patients with GD1 were included in this study. Of these, 308 had a heterozygous genotype in the GBA1 gene. Specifically, 113 patients carried the c.[1226A>G]; [1448T>C] genotype, while 195 had a heterozygous c.[1226A>G] genotype with a second variant distinct from c.[1448T>C]. Demographic data and clinical characteristics are detailed in Table 1.

When we analyze the group of patients with heterozygous c.[1226A>G] genotype with a second variant distinct from c.[1448T>C], the gender ratio is 1:0.8 (♂:♀), and almost half of the patients (55.41%) were diagnosed after the year 2000. The mean age at diagnosis was 24.7 years old (range: 1.0–77.0), and the mean age at start of therapy was 27.9 years old (range: 1.0–74.0). Thirty-eight patients (19.5%) were previously fully splenectomized.

During follow-up, 8 patients (4.8%) developed monoclonal gammopathy of undetermined significance (MGUS), 10 patients (5.1%) developed other malignancies, and 12 patients (6.1%) developed Parkinson’s disease. Of the 180 patients questioned, 20 patients (10.2%) were diagnosed with cholelithiasis, while 31 patients (17.2%) presented with other comorbidities. Fourteen patients (8.0%) had died by the time of the study. Regarding bone disease, 44 patients (22.51%) had advanced bone disease with skeletal complications, and 62 patients (50.8%) exhibited varying degrees of bone mineral density (BMD) loss.

In the Spanish Registry of GD, the most frequent allelic variant is c.[1226A>G] (75% of alleles), followed by c.[1448T>C] (18.5%). In Spanish GD patients, the most frequent genotype is heterozygous c.[1226A>G] (N370S) combined with c.[1448T>C] (L444P). However, 195 patients (52.13%) had a GBA1 heterozygous genotype other than NM_000157: c.[1226A>G] (N370S); [1448T>C] (L444P). These alternative genotypes are detailed in Table 2, with the new and classic protein nomenclature, and stratified according to functional severity.

2.2. Bone Involvement

A strong correlation between prior splenectomy and the presence of bone disease was observed (p = 0.005). In terms of genotype–phenotype correlation, patients with GBA1 genotypes other than c.[1226A>G]; [1448T>C] showed a significantly higher frequency of severe bone disease at diagnosis (p = 0.05) in terms as greater bone marrow infiltration and skeletal complications compared to patients homozygous for c.[1226A>G] and heterozygous for c.[1226A>G]; [1448T>C] (p = 0.017).

In our cohort, patients heterozygous for c.1226A>G with another variant (excluding c.1448T>C), as well as those with other genotypes (group B), exhibited significantly greater bone marrow infiltration and complications compared to patients with homozygous c.1226A>G (N370S/N370S) and heterozygous c.[1226A>G]; [1448T>C] (N370S/L444P) genotypes (group A) (p = 0.017) (Figure 1).

The N370S; N370S cohort included 53 patients, only 4 of whom were splenectomized. A few patients developed bone complications, and none developed PD.

The other group included 17 patients with generally mild variants with non-neurological symptoms, including mutations such as R463H, R496H, N188K, S13L, and R170C in a homozygous or heterozygous manner.

In addition, the severity of bone disease regarding gender was worse in males heterozygous for c.[1226A>G] with another variant (mean S-MRI: 10) vs. females (mean S-MRI: 6) (p < 0.001). In addition, there is a very clear relationship between bone mineral loss and S-MRI (p-value = 0.0001).

According to genotype, the presence of bone complications, such as infarct, necrosis, or fractures, was 59.9% for c.[1226A>G]; other different c.[1448T>C] vs. 38.3% for homozygous c.[1226A>G] and 41.0% for c.[1226A>G]; [1448T>C] (p = 0.002). These patients also had a significantly higher incidence of new bone crises during therapy (p = 0.0001). These groups demonstrated the highest degree of skeletal involvement. Moreover, patients homozygous for c.[1226A>G] (N370S) had a significantly lower prevalence of BMD loss compared to all other genetic subgroups (p = 0.003) (Figure 2). Moreover, a significant relationship was observed between bone mineral density loss and splenectomy patients (p = 0.0045), with a higher incidence in females (p = 0.023).

The predictive capacity of machine learning models to estimate the risk of severe bone involvement in Gaucher disease was further explored using a decision tree analysis (Figure 3). The model identified the S-MRI score as the most influential variable in determining bone status at diagnosis, with a threshold of 6.5 effectively separating patients with a higher probability of developing osteopenia or osteoporosis. Among those with lower S-MRI scores (≤6.5), age at diagnosis and genotype (excluding the common c.1226A>G or c.[1226A>G]; [1448T>C] variants) were key discriminators. In contrast, in patients with higher S-MRI scores, gender and age also contributed significantly to classification. The model’s discriminatory power supported the weight of the S-MRI score as the primary classifier (Table 3).

2.3. Development of Parkinson’s Disease

In the case of Parkinson’s disease, a significant relationship between the development of PD and genotypes other than homozygous c.[1226A>G] or heterozygous c.[1226A>G]; [1448T>C] is observed (p = 0.025). The very significant incidence of PD in patients with a family history of PD is noteworthy, independent of genotype (p = 0.00014). A relationship between PD and elevated ferritin levels can also be discerned (p = 0.02).

Despite the limited number of patients who developed Parkinson’s disease (n = 14), the decision tree analysis (Figure 4) yielded meaningful insights into the variables most predictive of this outcome. Age at diagnosis was the primary node, with younger age (≤44.5 years) being associated with lower risk. Among these patients, the absence of a family history of Parkinson’s disease and younger diagnostic age (≤3.5 years) further stratified those at reduced risk. Conversely, in patients diagnosed after 44.5 years of age, genotype emerged as a key discriminator: individuals with genotypes other than homozygous c.1226A>G or compound heterozygous c.[1226A>G]; [1448T>C] exhibited a substantially higher probability of developing Parkinson’s disease. Values across the tree support the relevance of genotype and family history as dominant variables (Table 4).

3. Discussion

Several international registries have been established to collect data from patients with Gaucher disease (GD), enabling the analysis of clinical, genetic, and therapeutic aspects across diverse populations. The first global registry (NCT00358943), launched in 1991 by the International Collaborative Gaucher Group (ICGG) with support from Genzyme-Sanofi, has provided valuable longitudinal data regardless of disease severity, treatment status, or therapeutic approach [18]. Similarly, the Gaucher Outcome Survey (GOS), funded by Shire-Takeda since 2010, has enrolled over 2000 patients as of early 2023, primarily from Israel and the United States, to assess the long-term safety and effectiveness of velaglucerase alfa and other GD therapies [19].

However, certain biases in the composition of these registries must be acknowledged. In GOS, nearly half of the patients were from Israel, where a higher prevalence of genotypes associated with mild disease phenotypes is well documented. A study by Zimran et al., which followed 103 untreated GD1 patients for over six decades, revealed a notably benign course, with mild hepatosplenomegaly, minimal rates of splenectomy, and no clinical bone involvement [20]. These findings underscore the significant influence of genotype on phenotypic expression and disease course.

In this context, Gaucher disease has emerged as a model within the lysosomal disorder group, owing to the availability of robust clinical and genetic datasets, validated biomarkers, and therapeutic diversity. Nowadays, the application of artificial intelligence (AI) and machine learning techniques to these datasets has proven especially valuable in identifying patterns and predicting clinical outcomes. Recent studies have successfully employed such methodologies in the diagnosis of rare diseases, in patient stratification using next-generation sequencing (NGS) [21], and in multi-omics integration [22].

Specifically in GD, machine learning models have shown high performance in identifying undiagnosed patients using real-world clinical and administrative data, significantly reducing diagnostic delays [23,24,25]. Our group has contributed to this field by integrating real-time registry data and biobank resources, allowing for in-depth analyses of clinical–biochemical correlations and the development of predictive models for disease-specific outcomes, such as skeletal complications [15,17], as well as for the identification of novel biomarkers [14].

Despite the extensive data generated, genotype–phenotype correlation in GD remains incomplete. This is partly due to the wide range of GBA1 variants, the differential genotype distribution across populations, and the possible epistatic effects of other genetic or epigenetic modifiers that remain poorly understood [26]. The Spanish population, with its rich genetic diversity and complex historical admixture, offers an opportunity to explore these associations [27]. In our cohort, the analysis of genotypes—particularly heterozygous combinations other than c.[1226A>G]; c.[1448T>C]—revealed a more severe clinical profile. This observation, not clearly established in previous international studies, highlights the importance of population-specific analyses.

Related to the discussion about other modifiers of the phenotype in GD, both non-modifiable alterations, such as other associated genetic alterations, and epigenetic modifiers are the subject of broad debate [26,28,29]. Some of the other concomitant external factors, such as infections, lifestyle habits, and environmental toxins, are modifiable, but these factors have not been included in this study. They could, however, be considered in future studies involving machine learning methods.

The application of AI-driven predictive models in our analysis further revealed genotype-associated clinical features, such as an increased risk of Parkinson’s disease and more severe skeletal involvement, including a higher prevalence of bone mineral loss, in certain genotypic subgroups [30,31]. Notably, male patients with these genotypes exhibited significantly greater skeletal burden compared to their female counterparts—a sex-specific trend that we previously reported, although its underlying pathophysiological basis remains unclear [17]. Sex-based differences in bone maturation have been determined via MRI images of healthy subjects, primarily in the lumbar spine [30]. Other relevant factors include the protective estrogenic effects on bone metabolism in females, together with gender-related differences in growth and bone remodeling [31]. Thus, the male sex may act as a modifier of skeletal severity in Gaucher disease. However, note that while this is evident in our cohort, where it is described for the first time, our conclusion remains speculative. These findings align with previous reports, such as that by Basiri et al. (2023) [18], which showed that patients with compound heterozygous GBA1 genotypes (e.g., c.[1226A>G]; other) had a higher risk of osteonecrosis compared to those with the homozygous c.[1226A>G] genotype [18].

Furthermore, the prevalence of osteopenia and osteoporosis was significantly higher among patients with c.[1226A>G]; other genotypes (p < 0.001), reinforcing the hypothesis of a more severe skeletal phenotype in this subgroup. These results are consistent with earlier observations suggesting that certain—albeit less common—genotypic combinations may confer a higher risk of bone disease, including osteonecrosis and chronic skeletal complications [29,31].

In contrast, no significant differences were observed in the incidence of malignancies across genotypic subgroups, suggesting that cancer risk in GD may be independent of GBA1 genotype [32].

However, our Random Forest analysis identified the c.[1226A>G]; other genotype as a strong predictor for the development of Parkinson’s disease, supporting growing evidence of genotype-associated neurological involvement. This is consistent with previous studies indicating an elevated risk of synucleinopathies in GD patients with compound heterozygous mutations, particularly those involving non-c.[1226A>G] alleles [33,34,35].

These insights underscore the value of combining detailed clinical registries with advanced computational approaches to unravel the multifactorial nature of Gaucher disease. Our data confirm that genotype plays a pivotal role in shaping clinical outcomes, particularly in the skeletal and neurological domains, and that male sex may act as an additional modifier of disease severity in specific genotypic contexts. Nevertheless, the scope and limitations of this study should be noted. All machine learning results only reflect internal validation on a single-center cohort (n = 195). We explicitly acknowledge that external generalization is limited and that estimates may still be optimistic despite cross-validation. Consequently, machine learning findings are interpreted as associative and hypothesis-generating; logistic regression remains the primary inferential analysis. Code and parameter settings used for training/validation are available from the corresponding authors upon reasonable request to facilitate reproducibility.

Finally, our findings highlight the importance of ongoing real-time data collection and multidisciplinary collaboration, including geneticists, clinicians, and computational biologists, in improving patient stratification and guiding personalized management strategies in rare diseases. Expanding this approach to larger and more diverse cohorts will be essential to validate these observations and further elucidate the underlying biological mechanisms.

4. Patients and Methods

4.1. Study Design and Population

The Spanish Gaucher Disease Registry (SGDR), established in 1996, collects demographic, clinical, genetic, laboratory, and imaging data on individuals diagnosed with Gaucher disease. To date, the registry includes data on 436 patients, of whom 374 have Type 1 Gaucher disease (GD1). For this study, we analyzed a subset of 358 GD1 patients with complete datasets.

All patients provided informed consent prior to inclusion in the registry. This study was conducted in accordance with Spanish biomedical research regulations (Law 14/2007) and the Declaration of Helsinki, and studies involving EU-based subjects must comply with the EU General Data Protection Regulation in the performance of all processing activities in connection with this study. All the patients included in the Spanish Registry have given their informed consent for the use of their data for research purposes. The scientific and ethics committees of Aragon (CEICA. PI11/00029 Date: 7 April 2011) and the FEETEG foundation approved this study. The blood samples were provided by the Biobank of the Aragon Health System, integrated in the Spanish National Biobanks Network (PT20/00112).

The anonymized dataset includes demographic data, clinical features at diagnosis, genetic findings, laboratory results, biomarker levels, and imaging assessments. Follow-up data include treatment regimens, occurrence of bone crises, development of malignancies and Parkinson’s disease, and survival outcomes.

4.2. Genetic Analysis

All patients underwent full GBA1 gene sequencing using Sanger sequencing to identify biallelic variants. In cases where one pathogenic variant was unidentified, long-range or nested PCR using gene-specific primers was performed to amplify the functional gene segments. These products were then sequenced to detect novel mutations. For patients with compound heterozygosity, the phase was determined by parental genotyping.

To rule out the 55-bp deletion in exon 9, all patients homozygous for the c.1226A>G variant were reanalyzed, as previously described [36]. Additionally, in patients apparently homozygous for c.1448T>C, two separate fragments (250 and 223 bp), encompassing exons 9 and 10, respectively, were sequenced to exclude pseudogene interference and confirm the presence of the variant within the functional gene.

4.3. Biomarker Analysis

Biomarker analyses were performed at the FEETEG reference laboratory. Chitotriosidase activity was measured using the method described by Hollak et al. (1994), with the artificial substrate 4-methylumbelliferyl-β-D-N,N′,N″-triacetylchitotrioside [37]. CCL18/PARC, an inflammatory chemokine produced by monocytes and T lymphocytes, was quantified using immunoassay techniques, as described by Boot et al. (2004) [38]. Glucosylsphingosine (GluSph), the deacylated form of glucosylceramide, which was first described in 1982 by Nilsson and Svennerholm [39], appearing to specifically reflect the accumulation of substrates in GD, was measured using a validated liquid chromatography–tandem mass spectrometry (LC-MS/MS) method [40].

4.4. Bone Involvement Assessment

Bone involvement was evaluated using S-MRI scoring [41] and a structured bone marrow report as an assessment tool for assessing bone marrow infiltration severity by MRI [42]. Evaluations were performed by a radiologist specialized in bone marrow skeletal imaging and stratified into three severity levels: mild, moderate, and severe. The bone mineral density was evaluated by DXA, and patients were classified according to the World Health Organization’s (WHO) criteria [43]. In addition, the overall clinical disease severity at diagnosis assessment was measured using the GD1 Severity Scoring System (GD1-DS3), including three major domains (i.e., bone, hematology, and visceral) [44].

4.5. Statistical Analyses

Descriptive statistics were used to summarize the data. Quantitative variables are reported as means with interquartile ranges, while categorical variables are expressed as percentages. Pearson’s correlation coefficient was used to assess linear correlations, and Chi-square tests were employed for categorical variables.

Data preprocessing, analysis, and modeling were conducted using R software (version 3.6.2), incorporating various packages such as car, ggplot2, vcd, GGally, plyr, igraph, rpart, dplyr [45,46,47]. For numerical explanatory variables, the Brunner–Munzel test was used, assessing associations by dividing the numerical data into two groups based on the binary variable categories. Fisher’s exact test was applied for binary variables, and the generalized Fisher’s exact test was applied for multi-categorical variables.

4.6. Predictive Modeling

Following the exploratory analysis, we used statistical and machine learning models to characterize associations between clinical variables and the target outcomes (bone mineral density loss and Parkinson’s disease). The primary inferential model was logistic regression (reporting odds ratios and 95% CIs). Decision tree-based models (single trees and Random Forest) were used secondarily to probe potential non-linearities and interactions that are not easily captured by linear terms [48,49]. Continuous predictors were kept in their original scale; no Winsorization was applied. Categorical variables were dummy-coded. We restricted model complexity to match sample size and outcome prevalence (events-per-variable principle), using penalized logistic regression (L2) when appropriate. For tree-based models, we constrained maximum depth, minimum samples per split/leaf, and the number of estimators to limit variance. Class imbalance (notably for PD) was handled with class-weighting in the loss function; any resampling was only performed within the training folds to avoid leakage. To enhance interpretability, we provide permutation-based variable importance for Random Forest, partial dependence profiles for the top predictors (S-MRI, age at diagnosis, genotype severity, sex), and simplified decision tree schematics that summarize the dominant splits in the Random-Forest ensemble. These visuals are intended to explain patterns rather than to represent the full model. More details are included in the Supplementary Materials.

5. Conclusions

Our study reinforces the pivotal role of the GBA1 genotype in shaping the clinical spectrum of Gaucher disease, particularly in relation to skeletal and neurological manifestations. Through the integration of real-world registry data, biobank resources, and machine learning models, we were able to identify specific genotypic combinations, especially those involving the c.[1226A>G]; other variants no c.[1448T>C], as strong predictors of disease severity and Parkinson’s disease risk.

These findings suggest the necessity of population-specific analyses, as genotypic distributions and their phenotypic consequences may vary significantly across cohorts. The Spanish population, characterized by high genetic diversity, provides a valuable context for uncovering less frequent but clinically relevant genotype–phenotype associations. Moreover, the observed sex-specific differences in skeletal involvement call for further investigation into potential biological modifiers of disease severity.

Overall, our results highlight the importance of combining clinical expertise with computational tools in the study of rare diseases. Ongoing data collection and interdisciplinary collaboration will be essential to refine patient stratification and advance personalized therapeutic strategies in Gaucher disease and other lysosomal disorders.

6. Highlights

Genotype–phenotype correlation remains incomplete in Gaucher disease, particularly for rare or compound heterozygous GBA1 variants.

Machine learning models identified c.[1226A>G]; other variants as strong predictors of Parkinson’s disease and severe skeletal involvement.

Sex-specific trends were observed, with male patients showing greater skeletal burden in certain genotypic subgroups.

Population-specific analyses revealed associations not previously reported in larger international registries.

Integrating clinical registries, biobanks, and AI tools enhances our ability to predict disease outcomes and tailor patient management.