Genetic modifiers of response to thalidomide in transfusion-dependent beta-thalassemia patients: a whole-exome sequence analysis

Introduction

Transfusion-dependent beta-thalassemia (TDT) patients require regular blood transfusions to maintain optimal hemoglobin levels (Musallam et al., 2021). The worldwide carrier rate of TDT is 1.5%. The highest carrier rate in the Southeast Asian region is found in Malaysia (12%), followed by India (8%) and then Pakistan (5–7%) (Khaliq, 2022; Shah et al., 2019). Stem cell transplant is the only cure, but a lack of compatible donors and its high cost are the problems (Cappellini et al., 2018; Musallam et al., 2021). Splenectomy, aimed at improving the red cell count and hence hemoglobin levels, is associated with sepsis and thrombosis, which makes it a less favorable option (Cappellini et al., 2018).

Current research in low-middle income countries (LMICs) is focused on discovering drugs that increase hemoglobin-F levels, rendering TDT patients transfusion-independent (Cappellini et al., 2014; Musallam et al., 2021; Taher et al., 2014). Thalidomide is an immunomodulatory drug that induces fetal hemoglobin (Musallam et al., 2021). However, the exact mechanism is yet to be discovered (Cappellini et al., 2018; Hagh et al., 2011). So far, researchers have only been able to report altered expression of genes in patients taking thalidomide (Hagh et al., 2011). Chen et al. (2017) has proposed that thalidomide increases levels of reactive oxygen species mediated by p38 MAPK signaling and histone H4 acetylation. However, this proposed mechanism fails to explain the raised hemoglobin-F by thalidomide. Yang et al. (2020) reported the association of certain single nucleotide polymorphisms (SNPs) in HBG2 and BCL11A gene with high fetal hemoglobin levels in thalassemia patients. Hagh et al. (2011) observed decreased expression of NF-kB in thalassemia patients taking thalidomide and hypothesized that thalidomide induces fetal hemoglobin by suppressing the expression of NF-kB. Liang et al. (2021) observed that the knockout of the GATAD2A gene caused increased levels of fetal hemoglobin and proposed that thalidomide acts by suppressing this gene. Despite extensive research, the exact mechanism of fetal hemoglobin induction by thalidomide remains unexplained.

Up to 30% of TDT patients do not respond to thalidomide (Bhattacharjee et al., 2023; Garg et al., 2023). Several possible mechanisms can exist. For instance, thalidomide pharmacokinetics and pharmacodynamics, the nature of the HBB mutation itself, or mechanisms related to gamma-chain repression and induction. Therefore, an unbiased discovery-based approach is more appropriate compared to a targeted approach. In this regard, whole-exome sequencing analysis would be a smart modality to precisely and comprehensively study all the genetic variations in protein-coding genes throughout the genome (Suhaimi et al., 2022).

In this first-of-its-kind study, we enrolled TDT patients and stratified them into excellent and non-responders to thalidomide. Genetic variations between the two groups were observed and documented as genetic modifiers or prognostic factors for response to thalidomide. They should be sequenced in all TDT patients before starting therapy to predetermine their response status.

Materials and Methods

This cross-sectional study was conducted as per the ethical principles outlined in the Declaration of Helsinki 1975 and after obtaining ethical approval from the Khyber Medical University (KMU) ethical committee (letter No: KMU/IPDM/IEC/2021/09). The study was conducted at the Institute of Pathology and Diagnostic Medicine (IPDM), Khyber Medical University (KMU), Peshawar from July 2021 to July 2023. Twenty TDT patients were recruited from two different centers i.e., “Blood Disease Clinic” and a “Thalassemia Clinic” using a non-probability purposive sampling technique. TDT patients of both genders and all ages who were given thalidomide for not less than 6 months’ duration were included while those with concomitant hematological disorders like leukemia, and those who had been transfused within 2 months from the date of determining their response status were excluded from the study. Compound heterozygous cases showing bands for hemoglobin-E, hemoglobin-C, and hemoglobin-S on electrophoresis were also excluded. The cases were followed up every month for up to 6 months to determine the response status. The laboratory investigations performed included complete blood counts, liver function tests, and renal function tests.

A written informed consent was obtained from the parents of the participants aged 1–6 years. Assent was obtained in cases where participants were aged 7–17 years of age. The research work was explained to the participants and their parents, ensuring them of the confidentiality of their data. They were informed that they could withdraw at any time from the research. They were also informed that the data will be disseminated for the purpose of adding knowledge to the existing literature and it will help in the treatment of thalassemia patients all over the world.

Results

A total of 20 TDT patients were included in the study. The mean age of study participants was 5.45 ± 3.81 years (range: 1–12 years). There were 17 (85%) males and 3 (15%) females, male to the male-to-female ratio being 5.6:1. Characteristics of study participants are given in Table 1.

In total, 222,180 germline variants were identified across 20 subjects, comprising 188,660 SNVs and 23,503 INDELs. To pinpoint functional mutations, stringent criteria were applied. Initially, splice site, nonsense, and frameshift mutations were isolated, followed by the inclusion of functional missense variants through bioinformatic annotation tools. From the refined variant pool, 24 candidate variants across 24 genes were identified (refer to File S5). Among these variants, 15 exhibited a nonsynonymous mutational burden of 0.99 in our 24 SNVs (Table 2).

The three most common SNPs in excellent responder group were found in exon region of CHI3L1 (rs880633), NPNT (rs35132891) and ZNF 208 (rs10425763) in 13/14 (92%) (p-value = 0.000), 12/14 (85%) (p-value = 0.001) and 10/14 (71%) (p-value = 0.011) excellent responders respectively. The commonest SNPs in non-responder group were found in exon region of PM20D1 gene (rs7518979) which was mutated in 5/6 (83.3%) (p-value = 0.000) non-responders, followed by LGR6 (rs75658797), MYH15 (rs4299484), and RESF1 (rs3207618) each of which were mutated in 4/6 (66.6%) (p-value = 0.003) non-responders (Table 2).

The functional enrichment analysis indicated that the genes of identified candidate variants in TDT were mainly enriched in biological processes related to thalassemia response such as IL-6, NF-kappa-B signaling, etc. (Fig. 1A). Identified rare variant-associated genes were subjected to pathway enrichment analysis based on gene ontology, KEGG (Fig. 1B). The majority of the pathways represented are associated with TDT such as Apoptosis, IL-17 signaling pathways, and B-cell receptor signaling pathways, lipid, folate, and cholesterol metabolism which are highly associated with TDT.

Discussion

To our knowledge, this is the first study that explored genetic variations associated with varying responses to thalidomide in TDT patients using whole-exome sequencing. In our study, genetic variations in NPNT, BMP3 & MAPK4 genes were seen in excellent responders. Currently, there is no established role of these genes in ineffective erythropoiesis in thalassemia. Additionally, in our study, genes known to participate in pro-apoptotic pathways and endoplasmic reticulum (ER) stress (TAF1C, INPP5D, ACP1) were mutated in the excellent responders, while the anti-apoptotic gene (BIRC2) was mutated in the non-responder group. This data hints towards genes that should be explored to understand the role of these genes as prognostic markers in the response of thalassemia patients to thalidomide.

The genetic variations in the genes highlighted in this study provide data for further transcriptome studies. This will help validate the association of these genes with response status to thalidomide in thalassemia patients and to explore the role of these genes in erythropoiesis.

The NPNT gene encodes an extracellular matrix protein, nephronectin, that binds to integrin α8β1. It regulates organogenesis and cell differentiation (Huang & Lee, 2005; Sopel et al., 2022). Recently, NPNT gene expression has been known to be regulated by the TNF-NFKB pathway and TGFβ-SMAD pathway (Fang et al., 2010; Sopel et al., 2022; Tsukasaki et al., 2012). However, the role of NPNT in erythropoiesis and hemoglobin synthesis is unknown (Magnussen et al., 2021). Clinical significance of SNP rs35132891 in the NPNT gene is not reported in Clinvar.

BMP3, the most abundant BMP in bone marrow, encodes a protein that binds to activin receptor type IIB, to activate SMAD2 which then regulates apoptosis and cell cycle (Wen et al., 2019). The clinical significance of SNP rs3733549 reported in our study is not reported in ClinVar. To our knowledge, there is no data available regarding the role of the BMP3 gene in erythropoiesis and thalassemia.

Disruption of the CHD-4 gene (components of the Nucleosome Remodeling and Deacetylase (NuRD) complex) is known to increase the fetal hemoglobin levels (Hossain et al., 2018). The SNP rs1639122 is reported as benign in ClinVar but in our study, it is associated with excellent response status. Only one study has shown a significant association of this SNP with bone mineral density (Liu et al., 2020).

The TAF1C gene encodes a protein that is involved in the synthesis of ribosomal RNA (rRNA), thus playing a role in processes like protein synthesis, cell differentiation, and apoptosis (Abdulwahab et al., 2019; Ferreira et al., 2020; Saproo, 2023). The role of the TAF1C gene in protein synthesis still needs further research (Saproo, 2023). We suggest research on the role of TAF1C in hemoglobin synthesis in thalassemia patients.

INPP5D encodes a 5′-phosphatase restricted to hematopoietic tissue. The role of INPP5D in erythropoiesis is unclear and needs further research (Mason et al., 2000). In our study, the SNP rs61752227 in the INPP5D gene was significantly associated with excellent response. However, the Clinical significance of SNP rs61752227 is not reported in ClinVar.

The SNP rs11553746 in the ACP1 gene as reported in our study, is reported as a missense variant in ClinVar. Recently, the role of the ACP1 gene in apoptosis and endoplasmic reticulum stress in HEPG2 cells has been reported (Bourebaba et al., 2020). Based on the finding that the ACP1 gene is mutated in excellent responders in our study, we suggest further studies to determine any possible role of the gene in ER stress-induced apoptosis in erythroid cells in thalassemia.

LGR6 receptors bind ligand R-spondins to activate the canonical Wnt-β-catenin pathway (Park et al., 2018). It is reported to be a tumor suppressor gene (Schuebel et al., 2007). Activation of LGR6 receptors on macrophages by maresin1 promotes phagocytosis and efferocytosis by macrophages (Chiang et al., 2019). Clinical significance of the SNP rs75658797 is not reported in ClinVar. However, we report that it is significantly mutated in non-responders.

CHI3L1 is found to be localized in the ER and its depletion, as seen in CHI3L1 knockout mice, increases superoxide dismutase-1 level which then activates ER-induced apoptosis in cells (Yu et al., 2023). We suggest further studies to evaluate any possible role of CHI3L1 in ER stress-induced apoptosis in ineffective erythropoiesis in thalassemia. In our study, the SNP rs880633 was present in excellent responders. However, its significance is not reported in ClinVar. However, the SNP rs880633 is reported to be associated with an increased risk of cystic fibrosis, chronic obstructive pulmonary disease, and hepatocellular carcinoma, while a decreased risk of developing deformities in patients with scoliosis (Hector et al., 2011).

ZNF208 is a member of the zinc finger gene family (Eichler et al., 1998). Generally, ZNF proteins bind to specific areas of DNA or RNA and regulate gene expression. We report the SNP rs10425763 in the ZNF208 gene to be associated with excellent response to thalidomide. However, the significance of the SNP is not reported in ClinVar. However, certain other variants of ZNF208 were reported to be associated with stroke and certain cancers (Alanazi & Iqbal, 2022; Wang et al., 2016; Yu et al., 2017).

The role of the MYH-15 gene in hematopoiesis is unknown and needs further workup. We report the SNP rs4299484 in the MYH15 gene to be associated with non-response status to thalidomide. A recent study established evidence for cross-talk between TGF-beta and MYH15 genes by showing that lung cancer cells, upon stimulation by TGF-beta showed increased expression of the MYH15 gene (Gladilin et al., 2019). This shows there is a link between the TGF-beta pathway and MYH15 gene expression.

RESF1 is an uncharacterized gene in humans. The gene is predicted to play a role in histone modification and gene silencing in association with the SETDB1 gene in mouse embryonic stem cells (Fukuda et al., 2018). The SNP rs3207618 is not reported in ClinVar. However, only one study reports this SNP to be associated with thalidomide-induced peripheral neuropathy in multiple myeloma patients (García-Sanz et al., 2017).

BIRC2 inhibits the non-canonical NF-KB pathway, a pathway that is necessary for the survival, proliferation, and differentiation of macrophages and other immune cells (Zadoroznyj & Dubrez, 2022). Loss of the BIRC2 gene is reported to lead to apoptotic cell death (Hrdinka & Yabal, 2019). Cellular stresses like hypoxia, ER stress, and DNA damage increase the expression of BIRC2 to promote apoptosis. However, research in erythroid cell lines in this regard has not been done.

The PM20D1 gene is involved in the regulation of body mass, temperature, insulin resistance, blood pressure, bone health, and angiogenesis. The precise functions of the body are however not fully known and need further research (Benson et al., 2019; Kim et al., 2020; Long et al., 2018). The gene has been reported as a quantitative trait locus associated with Alzheimer’s disease (Sanchez-Mut et al., 2018). There is no research regarding the function of the PM20D1 gene in hemoglobinopathies. We suggest studies to define the role of the P20D1 gene in thalassemia patients.

The major limitation of this study was the limited number of patients. This is attributed to the high cost of exome sequencing analysis. Baseline fetal hemoglobin levels were not monitored and compared in the two groups studied. A critical limitation of the study is the functional assay to confirm how these variants causally affect response to thalidomide. Further, inclusion of the pediatric population limits the generalization of findings to this subgroup only and not to the adult population.

We recommend larger studies to validate the findings of the study. Furthermore, functional studies like gene expression analysis or cellular assays should be done to establish the causal association of these variants with the response status to thalidomide.

Conclusions

This is the first study that shows that SNP rs880633, rs35132891, and rs10425763 are significantly associated with excellent response status to thalidomide in thalassemia patients, while rs7518979, rs75658797, rs4299484, and rs3207618 are associated with non-response status. We hypothesize for the first time that these genetic variants are the genetic modifiers of response to thalidomide in thalassemia patients and should be sequenced in all transfusion-dependent β-thalassemia patients before starting thalidomide to predetermine their response status. Only excellent responders will be continued on low-dose thalidomide therapy. However, gene expression analysis studies are recommended to establish these SNPs as the prognostic factors of response to thalidomide in thalassemia patients.

Supplemental Information