Genetics of glucocorticoid-associated osteonecrosis in children with acute lymphoblastic leukemia

Cure rates for childhood acute lymphoblastic leukemia (ALL) have approached 90% with therapeutic advances over the last several decades.^1-6  Intensification of therapy with glucocorticoids has played a crucial role in achieving these outcomes. However, one of the most common therapy-related and dose-limiting toxicities of therapy in children with ALL is glucocorticoid-induced osteonecrosis, particularly in those older than 10 years of age. The majority of symptomatic cases of osteonecrosis occur within the first 2 years of treatment,^7,8  often precipitating early withdrawal of glucocorticoids from therapy for ALL.

The incidence of glucocorticoid-induced osteonecrosis varies widely.^7,9  Age remains the most significant risk factor, with symptomatic osteonecrosis (defined as grades 2-4) occurring in 10% to 30% of children older than 10 years of age.^7,8,10-12  Glucocorticoid-induced osteonecrosis also complicates treatment of nonmalignant conditions such as solid organ transplant and arthritis.^13-15  Osteonecrosis can result in debilitation and adversely affect quality of life, often requiring surgical intervention.

In this study, we conducted the largest genome-wide association study (GWAS) to date of glucocorticoid-induced osteonecrosis, with replication cohorts including not only children treated for ALL⁷  but also adults and children treated with glucocorticoids for other medical conditions. Our goal was to identify germline genetic variants that predispose to glucocorticoid-induced osteonecrosis.

The discovery cohort included children with newly diagnosed ALL with germline DNA available who were treated on the Children’s Oncology Group (COG) AALL0232 protocol (NCT00075725) for high-risk B-precursor ALL (n = 2285) (Table 1; see supplemental Figure 1 and supplemental Table 1 on the Blood Web site). Validation cohorts included children with newly diagnosed ALL treated on the St. Jude (SJ) Total XV protocol (NCT00137111) (n = 361)⁷  (supplemental Figure 1 and supplemental Table 2), and a separate cohort comprising children and adults treated with corticosteroids in the Vanderbilt University Medical Center Biorepository BioVU database¹⁶  (n = 309) (Table 1; supplemental Figure 1 and supplemental Table 3). Patients included in the genetic association analyses represented 80% (n = 2285 of 2868) of participants on the COG AALL0232 protocol, and 73% (n = 361 of 498) of patients on the SJ Total XV protocol (supplemental Figure 1).

Informed consent was obtained from patients 18 years and older, and from parents or guardians for patients under the age of 18 in accordance with the Declaration of Helsinki. The COG AALL0232 protocol was approved by the National Cancer Institute and the Institutional Review Boards (IRBs) of participating institutions. The SJ Total XV trial was approved by the SJ IRB.

The second validation cohort for this study was derived from BioVU, Vanderbilt’s repository linking DNA from remnant blood samples to de-identified electronic medical record data.¹⁷  Due to de-identification, use of this resource qualifies as nonhuman subjects research. IRB exemption was required by the institution and was obtained prior to the study. Cases of glucocorticoid-induced osteonecrosis were identified using the de-identified electronic medical records in BioVU using keyword searches, International Classification of Disease, 9th edition codes, and exclusion criteria described in supplemental Materials. Controls were selected in a 3:1 (controls:cases) ratio, and matched on age, race, gender, and primary diagnoses. All cases and controls were manually reviewed to confirm inclusion/exclusion criteria (supplemental Figure 2).

Osteonecrosis was graded using the National Cancer Institute Common Terminology Criteria for Adverse Events, version 3.0 for SJ Total XV and version 4.0 for COG AALL0232, and categorized as absent (grade 0), asymptomatic (grade 1), moderate (grade 2), severe (grade 3), or disabling (grade 4). Patients with symptomatic grades 2-4 were considered to have osteonecrosis. For patients on AALL0232, patients with symptoms suggestive of osteonecrosis were evaluated by magnetic resonance imaging (MRI) to verify the diagnosis. In the SJ cohort, all patients were prospectively screened for osteonecrosis with serial MRI, regardless of symptoms. Details for case ascertainment in the Vanderbilt cohort are provided (supplemental Materials and methods).

Single nucleotide polymorphism (SNP) genotyping was performed on germline DNA from patients from the COG AALL0232 and SJ Total XV protocols using the Illumina Human Exome BeadChip v1.1 and Affymetrix GeneChip Human Mapping Array 6.0. The Illumina Exome Chip, Omni1-Quad, 1M-Duo, and Human660W-Quad BeadChip Arrays were used for the Vanderbilt BioVU cohort. For SNPs not interrogated on the arrays, imputation was performed using the 1000 Genomes Project (http://www.1000genomes.org/) database as the reference genome, leveraging linkage disequilibrium (LD) within racial ethnic groups (Northern European ancestry, West African ancestry, and other) with MaCH-Admix software (University of North Carolina) for the COG AALL0232 and SJ Total XV cohorts, and IMPUTE2 (University of Oxford, United Kingdom) for the Vanderbilt cohort.^18-20  For the imputation, reference groups for imputing SNPs in patients classified as white (see below) were as follows: the European individuals in 1000 Genomes, for patients classified as black were the 1000 Genomes Africans, and the remaining patients were all individuals in the 1000 Genomes database.

Genetic ancestry for each patient was determined using STRUCTURE (version 2.2.3).²¹  When ancestry was assessed as a categorical variable, individuals were classified as white, black, Hispanic, and Asian based on percentage-inferred genetic ancestry as follows: >90% Northern European (CEU) were classified as white; >70% West African (YRI) were classified as black; >90% East Asian (CHB/JPT) were classified as Asian; those with Native American ancestry >10% and greater than the percentage of West African ancestry were classified as Hispanic. Patients not falling into these groups were categorized as “Other.”

Identical quality control measures were enforced for the discovery and validation cohorts. All SNPs with a call rate <95% were excluded. SNPs with a minor allele frequency (MAF) >1% but deviating from Hardy–Weinberg equilibrium within Europeans (P < .0001) were also excluded. Otherwise, no MAF threshold was enforced.

For the discovery GWAS, SNP genotypes were compared in 250 ALL osteonecrosis cases and 2035 controls enrolled on COG AALL0232. Adjusting for gender, age, percent ancestry as a continuous variable, and treatment (see supplemental Materials and methods for details on consideration of treatment variables), association of genotypes with ON was tested with a Cox proportional hazard model for time-dependent analyses and logistic regression for time-independent analyses. For the time-independent analyses, only patients with a follow-up time of 800 days or greater from the start of therapy on COG AALL0232 were included in the analysis. Results from analyses with imputed SNPs and each independent platform were merged and rank ordered by P value. Analyses were performed using R software (version 3.0; www.r-project.org).

We excluded rare/low frequency SNPs (MAF <0.1) with a protective, negative correlation with osteonecrosis risk (odds ratio [OR] <1).

Time-independent analyses were performed using a logistic regression model in the 2 validation cohorts: 68 osteonecrosis cases and 293 controls from the SJ Total XV cohort as previously published,⁷  and 82 cases and 227 controls from the Vanderbilt BioVU cohort. Meta-analysis of the discovery and validation cohorts was performed using Stouffer’s Z-score method.²² 

Pathway and network analyses were performed to identify biological networks enriched within the top ranked SNPs for each of the cohorts using QIAGEN Ingenuity Pathway Analysis (www.qiagen.com/ingenuity). Top SNPs were selected with meta-analysis P value cutoffs < .0001 or < .001, and genes closest to SNPs based on chromosomal location were used for the pathway analyses.

The phenome-wide association studies (PheWAS) analysis was performed as described^23,24  (supplemental Materials and methods) for the 137 SNPs within the glutamate receptor signaling pathway genes. We tested for their associations with 1358 phenotypes using the extant BioVU genotype data linked to electronic health records for over 10 000 individuals.

The overall frequency of symptomatic grades 2-4 osteonecrosis in patients included in the entire COG AALL0232 cohort (all ages) was 10.9%. For the discovery COG AALL0232 cohort, a Cox proportional hazard model with nongenetic and clinical factors was used to identify covariates to include in the GWAS. Multivariate analysis revealed age ≥10 years (hazard ratio [HR] = 11.1; 95% CI, 6.7-18.1; P = 1.57 × 10⁻²¹), and female gender (HR = 1.39; 95% CI, 1.08-1.78; P = 9.82 × 10⁻³) (supplemental Table 1) to be associated with a higher risk of osteonecrosis. African genetic ancestry was associated with a decreased risk of osteonecrosis compared with individuals of European genetic ancestry (HR = 0.19; 95% CI, 0.07-0.51; P = 1.02 × 10⁻³). Treatment variables A and B (supplemental Materials and methods and supplemental Table 4) were also included as covariates in the GWAS.

As reported previously,⁷  the cumulative incidence of symptomatic osteonecrosis grades 2-4 in the SJ Total XV cohort was 17.6%, with older age (P = 7.81 × 10⁻⁷) and more intensive therapy (P = 9.48 × 10⁻³) as significant covariates (supplemental Table 2).⁷ 

In the discovery GWAS in the COG AALL0232 cohort, adjusting for gender, age, ancestry, and treatment, a variant at 9q31.1 (rs10989692, P = 3.59 × 10⁻⁷), ∼170kb 5′ of the GRIN3A locus (glutamate [NMDA] receptor subunit 3A), had the strongest association with osteonecrosis risk in both a time-dependent and time-independent analysis (supplemental Tables 5 and 6; supplemental Figure 3). Patients with an additional A allele at rs10989692 had a higher risk of osteonecrosis (HR = 2.03; 95% CI, 1.55-2.66) (Figure 1A), with 73 of 250 cases carrying a risk allele (69 heterozygous and 4 homozygous). The correlation between rs10989692 genotypes and osteonecrosis was stronger in patients with >90% European ancestry (n = 1268; HR = 2.22; 95% CI, 1.64-2.99) (Figure 1B), risk allele frequency (RAF) in patients with and without osteonecrosis (0.173 vs 0.089; P = 1.7 × 10⁻⁶) than in patients of non-European ancestry (n = 1017; HR = 1.63; 95% CI, 0.98-2.63; P = .057) (Figure 2; supplemental Figure 4). SNP rs10989692 was further validated in the Vanderbilt cohort (OR = 2.26; P = .0074), again with a stronger association in whites (n = 216; OR = 2.68; P = .012; RAF: 0.132 in cases vs 0.06 in controls) than in those with non-European ancestries (n = 93; OR = 1.72; P = .26). The effect of rs10989692 in the SJ cohort was marginal when all patients were included (OR = 1.87; P = .063), but reached statistical significance when only patients with >90% European ancestry were considered (n = 260; OR = 2.29; P = .028; RAF: 0.123 in cases vs 0.086 in controls). The RAF was very similar in both the discovery and replication cohorts (0.106, 0.116, and 0.130 for COG, SJ, and Vanderbilt cohorts, respectively). An additional SNP (rs28584318) in high LD with rs10989692 (r² = 1 and 0.83 in HapMap CEU and African populations, respectively) was also associated with osteonecrosis (P = 4.86 × 10⁻⁷; P = .07; P = .0073 in AALL0232, SJ, and Vanderbilt cohorts, respectively).²⁵  We evaluated the effect of the rs10989692 SNP across different ancestral groups by analysis of variance test of the interaction term of the genotype and genetic ancestry in AALL0232. The P value for the interaction term, P = .075, suggests that there may be a difference in the effect size of this polymorphism across ancestral groups. When restricting the analysis to patients >10 years, the risk associated with rs10989692 was maintained in AALL0232 (n = 1468; HR = 2.07; 95% CI, 1.59-2.70; P = 1.44 × 10⁻⁶; RAF 0.152 in cases vs 0.0905 in controls) with a similar trend in the SJ cohort that did not reach statistical significance (n = 91; HR = 1.67; 95% CI, 0.55-5.10; P = .37; RAF 0.10 in cases vs 0.077 in controls) (supplemental Tables 7 and 8; supplemental Figure 3C), and interaction analyses between genotype and age did not show a significant interaction between GRIN3A rs10989692 genotype and age (P = .86 and P = .89) in the AALL0232 and TOTALXV cohorts, respectively.

Meta-analysis was performed by combining GWAS results from both the discovery cohort and the 2 validation cohorts. The 9q31.1 locus (near GRIN3A) was the highest ranked region, including multiple SNPs in LD with rs10989692 (P = 2.68 × 10⁻⁸) (Table 2; Figure 3). After the 9q31.1 locus, the next highest ranked variant was an intronic SNP rs2154490 in glutamate receptor, ionotropic, kainate (GRIK1) on chromosome 21 (P = 1.28 × 10⁻⁶) (Table 2; Figure 3). The additional minor A allele at rs2154490 conferred a higher risk of osteonecrosis in all 3 cohorts: COG AALL0232 (HR = 1.29; P = .016), SJ Total XV (OR = 1.86; P = .0078), and Vanderbilt (OR = 2.13; P = 9.1 × 10⁻⁴) (supplemental Figure 5).

The top nonsense SNP from the meta-analysis (supplemental Table 9) was within pleckstrin homology domain containing, family H (with MyTH4 domain) member 1 (or PLEKHH1) (P = .0023), involved in phosphate binding and implicated in liver lipid homeostasis.²⁶  Another nonsense SNP (P = .0056) was within STEAP4, a ferrireductase involved in osteoclastogenesis²⁷  and adipocyte development.²⁸  Additional missense loci of interest from the meta-analysis include those in ZFHX3, a gene involved in atrial fibrillation and Kawasaki disease²⁹  (P = 1.46 × 10⁻⁶), and loci in COL22A1 (collagen, type XXII, α 1; P = .0003), involved in maintaining integrity of tissue junctions.

Based on the meta-analysis, there were 197 SNPs with P values < .0001, which were annotated to 64 genes (supplemental Table 10). Ingenuity Pathway Analysis using these 64 genes (supplemental Table 11) showed that the glutamate receptor-signaling pathway, including 3 genes GRIN3A, GRIK1, and GRM7, was the top canonical pathway (P = 4.8 × 10⁻⁴). Using a P value cutoff of .001, there were 433 genes, and the glutamate receptor-signaling pathway remained the top pathway, with additional genes in the pathway including GRM3, GRIK4, and GRIA1 (P = 9.8 × 10⁻⁴) (supplemental Table 11).

All SNPs with P value < .05 in or near genes within the glutamate pathway found to be significant in at least 1 of 3 cohorts were selected for analysis in a PheWAS (supplemental Materials) of the Vanderbilt University Medical Center BioVU database to explore whether these same osteonecrosis variants were associated with any of 1358 other phenotypes in the BioVU repository (supplemental Table 12). None of these associations achieved a Bonferroni significance threshold of 3.7 × 10⁻⁵. However, among the top 5 phenotypes associated with SNPs in GRIN3A were diseases of the respiratory tract (OR = 3.02; P = 3.1 × 10⁻⁴) and the long-term use of antithrombotics within the circulatory system (OR = 2.13; P = 4.7 × 10⁻⁴) (supplemental Table 13). Other phenotypes included cerebral ischemia (OR = 1.64; P = 2.5 × 10⁻³), and arterial embolism and thrombosis (OR = 1.88; P = 4.2 × 10⁻³).

The development of osteonecrosis can result in serious debilitation and requires surgical intervention in many cases. For example, in the Children’s Cancer Group 1961 study, 143 patients (of 2056 enrolled) developed symptomatic osteonecrosis at 377 confirmed sites leading to 139 surgeries.⁸  Identification of risk factors for the development of osteonecrosis in patients undergoing therapy for ALL might facilitate tailoring therapy to minimize its risk.

In this study, we performed a GWAS, which identified a locus at 9q31.1 near GRIN3A as associated with osteonecrosis in both the discovery and validation cohorts, and the second highest ranked variant was in a related gene, GRIK1. The same GRIN3A variant was associated in the meta-analysis of the discovery and validation cohorts (P = 2.68 × 10⁻⁸), and the glutamate receptor pathway was the top canonical pathway based on the meta-analysis of all 3 cohorts. These findings suggest the involvement of the glutamate pathway in the pathogenesis of osteonecrosis not only in childhood ALL, but in a heterogeneous cohort of children and adults who were treated with prolonged corticosteroids for other medical conditions.^30,31  Interestingly, we previously reported that polymorphisms within a third glutamate receptor gene, GRIA1, were the top genomic variants associated with asparaginase allergy in a GWAS of children with ALL in the SJ cohort,³²  a finding recently independently replicated.³³  We have also reported that asparaginase allergy was associated with lower systemic exposure to asparaginase,³⁴  which in turn was associated with a lower risk of osteonecrosis in the SJ cohort.³⁵  This prompted us to evaluate asparaginase allergy in the COG AALL0232 discovery cohort (n = 1845 with follow-up >800 days), and consistent with the SJ cohort, grades 2-4 asparaginase allergy was associated with the risk of osteonecrosis: 206 out of 1647 (12.5%) for those without asparaginase allergies had osteonecrosis compared with 12 out of 198 (6.7%) for those with asparaginase allergies (P = .0096).

Concomitant asparaginase treatment is associated with higher plasma exposure to dexamethasone.^7,34  Although it is possible that the glutamate receptor pathway association was via secondary effects on asparaginase, which affects risk of osteonecrosis in ALL patients,³⁵  this does not explain the association of the glutamate pathway within the BioVU cohort, whose patients did not receive asparaginase, and thus also implicates a mechanism with a direct effect of glutamate receptors. Multiple mechanisms can lead to osteonecrosis, with contributions to risk from therapy and host-specific factors.

Previously described biologic mechanisms by which glucocorticoids may induce osteonecrosis include thrombosis, hyperlipidemia-associated enlarged adipocytes in bone, arteriopathy, and direct toxicity to osteocytes.^36-41  Prior candidate gene and genome-wide investigations have implicated polymorphisms involved in lipid homeostasis (ACP1), fibrinolysis (SERPINE1), antifolate pharmacodynamics (thymidylate synthetase), and glucocorticoid response (VDR).^7,10,42,43  These variants did not replicate in the discovery cohort, possibly due to differences in therapy and methods of case/control ascertainment. Differences in therapy and in frequency of performing MRIs may also explain the different rates of osteonecrosis seen in the AALL0232 and SJ Total XV protocols. Importantly, the GRIN3A rs10989692 variant was associated with osteonecrosis in both cohorts despite these underlying protocol differences.

Glutamate receptor genetic variations have not been previously reported as risk factors for osteonecrosis. Glutamate is released by osteocytes in response to mechanical load, which opens stretch-sensitive calcium channels and causes activation of osteoblast receptors,⁴⁴  and glutamate impairs endothelial barrier function.⁴⁵  Glucocorticoids have been shown to induce the expression of glutamine synthetase in osteoblasts⁴⁶  and hepatoma⁴⁷  cells. Genetic variation in GRIN3A has been associated with the severity of vascular complications of Kawasaki disease.⁴⁸  Interestingly, the disruption of the vascular supply to bone is a proximal event to glucocorticoid-induced osteonecrosis in a murine model, with or without asparaginase.⁴¹  Herein, in our PheWAS of GRIN3A SNPs, vascular phenotypes including cerebral ischemia, arterial embolism, and thrombosis trended toward association with glutamate receptor variants. Thus, glutamate and variations in glutamate receptors may contribute to a proximal vascular event that leads to an increased risk of osteonecrosis in individuals exposed to glucocorticoid therapy.

In conclusion, we hypothesize that different mechanisms of glucocorticoid-induced osteonecrosis may predominate among patients, influenced by variation in concurrent drug therapy as well as inherited genetic risk factors. To our knowledge, this is the largest genome-wide investigation of glucocorticoid-induced osteonecrosis. Our findings suggest for the first time a possible association between inherited genetic variations in glutamate receptors and the development of glucocorticoid-induced osteonecrosis in both ALL and non-ALL settings. The ability to identify genetic risk factors for osteonecrosis has implications for understanding the underlying mechanism of this common and serious adverse effect of glucocorticoids, and may have implications for modifying therapy decisions in the future.

This study was supported by the National Institutes of Health, Institute of General Medical Sciences, National Cancer Institute, Human Genome Research Institute, and National Library of Medicine grants GM 92666, GM 115279, CA142665, CA 21765, CA 36401, CA98543 (COG Chair’s grant), CA98413 (COG Statistical Center), CA114766 (COG Specimen Banking), U01-HG04603, RC2-GM092618, and R01-LM010685; the Vanderbilt Clinical Translational Science Award grants UL1-TR000445, 5T32-GM007569, and LLS 6168; and the American Lebanese Syrian Associated Charities.

Contribution: E.A.R., C.-H.P., and M.V.R. were responsible for concept and design; S.L.V.D., S.K., E.B., D.M.R., J.C.D., E.L., N.W., M.L.L., S.P.H., P.S., M.D., L.A.M., and M.V.R. collected and assembled the data; W.Y., S.E.K., T.Y.C., M.B., L.B., D.M.R., J.C.D., W.L.C., C.C., D.P., C.A.F., C.L., C.S., P.S., S.J., W.E.E., M.D., L.A.M., and M.V.R. were responsible for the analysis and interpretation of data; and all authors contributed to the writing of the manuscript.

Conflict-of-interest disclosure: The authors declare no competing financial interests.

Correspondence: Mary V. Relling, Pharmaceutical Department, St. Jude Children’s Research Hospital, 262 Danny Thomas Pl, Memphis, TN 38105; e-mail: mary.relling@stjude.org.