Genetic risk factors for the development of osteonecrosis in children under age 10 treated for acute lymphoblastic leukemia

Progressive intensification of multiagent chemotherapy has improved outcomes for children with acute lymphoblastic leukemia (ALL), with cure rates exceeding 90%.^1-6  Results have been particularly favorable for patients <10 years of age with National Cancer Institute (NCI) standard-risk (SR) ALL.^7,8  Unfortunately, such intensified therapy has also been associated with significant increases in the occurrence of therapy-related osteonecrosis.^9-12 

Prior studies have identified multiple clinical risk factors for the development of osteonecrosis, including female sex, European ancestry, the administration of 3 weeks of continuous rather than alternate-week dexamethasone during delayed intensification, and intensive vs standard therapy.^9-11  However, age remains the strongest and most consistently identified factor, with patients 10 to 20 years old at greatest risk.^10,11,13-15 

Given this age-related risk, the majority of children in prior investigations of genetic predisposition to osteonecrosis were >10 years.^10,15  However, because ALL is so common in young children, up to 40% of osteonecrosis cases develop in children <10 years of age.¹⁰ 

In this study, we identified genetic factors associated with the development of osteonecrosis in the largest cohort of SR ALL patients evaluated to date and validated these findings in a cohort of NCI high-risk ALL patients <10 years of age.

The discovery cohort consisted of children with newly diagnosed SR B-precursor ALL treated on the Children’s Oncology Group (COG) AALL0331 (NCT00103285) protocol with germline DNA available; of these, all patients with osteonecrosis as of June 30, 2012 were submitted for genotyping as cases. Controls were taken from a previously genotyped subcohort of AALL0331 (supplemental Methods, available on the Blood Web site). Of 111 cases who were evaluable after induction and developed grade 2 to 4 osteonecrosis, genotyping was completed for 96. Patients (N = 34, 14 cases and 20 controls) who received alternate-week dexamethasone after protocol amendment 2C (see supplemental Methods for details) were excluded from the analysis (Figure 1). Controls were excluded if <1000 days of follow-up data were available, a time when the majority of osteonecrosis was detected. Following these steps, the discovery cohort included 82 osteonecrosis cases and 287 controls (Figure 1).

The replication cohort consisted of children <10 years of age at diagnosis treated for newly diagnosed high-risk B-precursor ALL on the COG AALL0232 (NCT00137111) protocol. The characteristics of this cohort have been previously described.¹⁵  The analysis was limited to the 817 children (20 cases, 797 controls) with genotype and phenotype data available (supplemental Figure 1).

Informed consent was obtained from patients ≥18 years of age and from parents or guardians of patients <18 years of age in accordance with the Declaration of Helsinki. The COG AALL0331 and AALL0232 protocols were approved by the NCI and the institutional review boards of participating institutions.

Patients on both protocols were prospectively monitored for clinical signs and symptoms of osteonecrosis, and the diagnosis was confirmed by imaging according to institutional preference. Osteonecrosis was considered a reportable event and was graded using the NCI Common Terminology Criteria for Adverse Events Version 4.0; controls had absent (grade 0) or asymptomatic (grade 1) osteonecrosis; cases had moderate (grade 2), severe (grade 3), or disabling (grade 4) osteonecrosis.

Genotyping of single nucleotide polymorphisms (SNPs) for AALL0331 and AALL0232 was performed using the Illumina Human Exome BeadChip v1.1 and the Affymetrix Gene Chip Human Mapping Array 6.0. For SNPs not interrogated by direct genotyping, genotypes were imputed using MaCH-Admix software (University of North Carolina), using phase 1 release v3 of the 1000 Genomes Project (http://www.1000genomes.org) as the reference genomes.¹⁶ 

For a subset of 15 cases of osteonecrosis, whole exome sequencing was performed. Library generation was performed using NimbleGen SeqCap EZ Exome Enrichment Kit v2.0, and sequencing was performed using Illumina HiSequation 2000. Raw reads were aligned against human genome hg19 using Burrow-Wheeler Aligner.¹⁷  Raw reads were converted to genotypes using the Genome Analysis Toolkit version 3.2.¹⁸ 

Genetic ancestry for all patients was defined using STRUCTURE v2.2.3 as previously described.¹⁹  When ancestry was assessed as a categorical variable, individuals were classified as white, black, or Hispanic based on percentage inferred genetic ancestry as follows: >90% Northern European (CEU) were classified as white; >70% West African (YRI) classified as black; those with Native American ancestry >10% and greater Native American than West African ancestry were classified as Hispanic. Patients not falling into these groups were categorized as Other.

All SNPs with a call rate of <95% or poor clustering were excluded. SNPs were excluded if they deviated from Hardy-Weinberg equilibrium in whites (P < 1 × 10⁻⁴). SNPs with a minor allele frequency (MAF) >5% and SNPs with a MAF >0.5% with the minor allele as the risk allele for the development of osteonecrosis were included.

AALL0331 treatment arms were clustered into therapy groups based on premaintenance therapy duration, dexamethasone and asparaginase dosing, and previously reported therapy groups at increased risk for the development of osteonecrosis (supplemental Methods; supplemental Table 1).²⁰  Within the discovery cohort, SNP genotypes were compared between cases and controls while controlling for age, genetic ancestry as continuous variables, sex, and therapy group. Evaluation used a time-to-event proportional hazard regression analysis with patients censored at the time of diagnosis of osteonecrosis or at last follow-up. Relapse and off-protocol therapy were treated as competing events. Analysis was performed using the R survival package (v2.38-1) implemented in PLINK (v1.07; http://pngu.mgh.harvard.edu/purcell/plink/) using Rserve.^21,22  Additional data analysis was performed using R (v3.1.1; www.r-project.org).²³ 

To test for validation, a parallel genome-wide association study (GWAS) was performed in the replication cohort. SNP alleles identified in the discovery cohort were considered validated if the P values in both cohorts were <.05 (supplemental Figure 2). A meta-analysis was performed for all SNPs passing quality control in both the discovery and replication cohorts using Stouffer’s Z-score method weighted by the GWAS standard error of the SNP in each cohort.²⁴ 

Enhancer enrichment analysis was performed using HaploReg (v2).²⁵  Validated SNPs with a meta-analysis P < .001 were evaluated for enhancer enrichment across available Roadmap samples. Associations were sorted according to their enrichment for enhancers and strongest enhancers.²⁶ 

SNPs validated in both cohorts with a meta-analysis P < .001 were annotated to genes, with the closest gene selected for intergenic SNPs. Pathway analysis was performed using QIAGEN Ingenuity Pathway Analysis (www.qiagen.com/ingenuity) to identify canonical biological pathways enriched within top SNPs. Only pathways with ≥5 genes associated with osteonecrosis were considered.

We identified covariates to include in the GWAS for osteonecrosis in the discovery cohort (supplemental Table 2), which included age (P = 1.88 × 10^-10), therapy group (P = 6.8 × 10⁻⁴ for the difference across all 4 groups), sex (P = .008), and genetic ancestry (P = .02 for European ancestry).

A GWAS of the discovery cohort adjusting for the covariates of age, sex, therapy group, and ancestry (supplemental Table 2) identified 16 SNPs in 6 genes associated with osteonecrosis at a P value threshold of <5 × 10⁻⁸ (supplemental Table 3). Two SNPs, rs77556622 and rs76599360 (both at P = 1.13 × 10⁻⁹), were located ∼11 and 20.7 kb, respectively, 3′ of BMP7 (bone morphogenic protein 7). The risk alleles (C and T, respectively) for both SNPs were present in 5 of 82 cases and 0 of the 287 controls and conferred a 22-fold increase in osteonecrosis risk (95% confidence interval [CI] 8.15-59.6). Patients with these SNPs developed osteonecrosis despite clinical characteristics placing them at low risk, including 3 patients in the lowest risk treatment group, 2 males, and 1 patient with black ancestry (Figure 2).

In the replication cohort, 10 of the top 100 SNPs from the discovery cohort were validated, including 2 SNPs (rs79085477, rs75161997) in BMP7 in full linkage disequilibrium (LD) with one another, with P = 5.34 × 10⁻⁸ (hazard ratio [HR] 15.0, 95% CI 5.64-39.7) in the discovery cohort and P = .0498 (HR 8.43, 95% CI 1.002-71.1) in the replication cohort. The top replicated SNPs represented 5 loci across 4 genes including BMP7 (2 SNPs), PROX1-AS1 (6 SNPs), LINC00251 (1 SNP), and DOK5 (1 SNP) (Table 1). Interestingly, the top validated SNPs in both BMP7 and PROX1-AS1 (rs7516997 and rs1891059, respectively) are predicted to increase glucocorticoid receptor binding based on binding site motif alterations.

Meta-analysis of both cohorts identified 5 loci with genome-wide significant associations with osteonecrosis (P < 5 × 10⁻⁸), 3 of which were significant in both cohorts, including loci near PROX1-AS1 and BMP7 (Table 1; Figure 3). The top SNP from the meta-analysis (rs9632544 near ERV3-1) was significant in the discovery but not the replication cohort (P = 2.36 × 10⁻⁹ and P = .052, respectively).

In the meta-analysis, the top replicated coding SNP was rs34144324, a missense variant in GRID2 (glutamate receptor, ionotropic, δ2) with a meta-analysis P = 7.07 × 10⁻⁷ (HR 3.78, 95% CI 2.23-6.39; P = 8.65 × 10⁻⁶ and P = .014 in the discovery and replication cohorts, respectively). In the 15 patients with whole exome sequencing available, 12 nonsynonymous variants were sequenced with a P < .001 in the meta-analysis. There was 99.4% agreement in variant calls between the array-based genotype calls and the sequencing genotype calls, including 100% agreement for the GRID2 SNP.

SNPs validated in both cohorts with a meta-analysis P < .001 (N = 3271) were evaluated for enhancer enrichment using HaploReg. We focused on tissues in which SNPs were enriched in both all enhancers and the strongest enhancers: colon smooth muscle, adipose-derived mesenchymal stem cells (MSCs), bone marrow–derived MSCs, and duodenal smooth muscle. When restricting the analysis to validated SNPs with a meta-analysis P < 1 × 10⁻⁵ (N = 140; supplemental Table 4), enhancer enrichment was observed only in bone marrow–derived MSCs (2-fold enrichment for all enhancers [P = .0148] and 4.3-fold enrichment for the strongest enhancers [P = 2.81 × 10⁻⁴]). Adipose-derived MSCs remained enriched for SNPs in all enhancers (1.8-fold, P = .011), with marginal enrichment in the strongest enhancers (1.9-fold, P = .068; Table 2).

There were 3271 SNPs with meta-analysis P < .001 and significant associations with osteonecrosis (P < .05) in both the AALL0331 and AALL0232 cohorts (the top 140 with P < 1 × 10⁻⁵ are shown in supplemental Table 4). The 459 genes to which these SNPs annotated were significantly enriched within 8 canonical pathways. Of these, 7 contained overlapping glutamate receptor genes, with the glutamate receptor signaling pathway most overrepresented (P = 9.92 × 10⁻⁴, 10.5% of genes in pathway associated with osteonecrosis). The adipogenesis pathway was the only overrepresented pathway (P = .0148, 5.5% of genes in pathway associated with osteonecrosis) whose genes did not overlap with glutamate receptor signaling (Table 3; supplemental Table 5).

Osteonecrosis has been recognized in recent years as a major toxicity of pediatric ALL therapy. Children <10 years of age constitute ∼75% of pediatric ALL cases.²⁷  Thus, although osteonecrosis incidence is low in this age group, those <10 constitute up to 40% of cases of osteonecrosis.¹⁰  It is hypothesized that inherited risk factors may be more penetrant in young rather than older patients. Better understanding of the individual risk factors for the development of osteonecrosis in younger children and the underlying biology in this population and in older patients is needed to improve prediction and management of this often debilitating complication. These patients who developed osteonecrosis despite being at low risk based on age and treatment factors may provide unique insights into the genetic variations that influence development of osteonecrosis.

The occurrence of osteonecrosis is highly dependent on the therapeutic regimen that patients receive. Thus, the identification of genetic variants that increase osteonecrosis risk in patients <10 years of age treated on both SR (AALL0331) and high-risk therapy (AALL0232) protocols suggests the identified variants may be important to osteonecrosis risk in this age group in the context of any modern ALL protocol (using both minimal residual disease–directed intensification of therapy, pegylated-asparaginase, and a delayed intensification phase).

We performed the first GWAS for osteonecrosis risk in patients exclusively <10 years of age and identified novel risk variants near BMP7 and PROX1-AS1 (Table 1). Studies have shown that BMP7 is released in response to bone damage in osteoarthritis and spondyloarthritis, and its release is increased by mechanical stress.^28-31  Bone morphogenic proteins can induce mesenchymal precursor cells to differentiate into osteoblasts and inhibit the formation of osteoclasts,^32,33  and thus variants affecting this gene could contribute to osteonecrosis by altering bone metabolism and formation both prior to and during ALL therapy. Moreover, BMP7 is known to be toxic to vascular smooth muscle.³⁴  The prominent role of local arteriopathy, at least in part due to vascular smooth muscle damage, in the development of osteonecrosis³⁵  supports the possibility that BMP7 could contribute to osteonecrosis via direct toxic effects on local bone vasculature.

PROX1 has been shown to control the differentiation of lymphatic endothelial cells from vascular endothelial cells.³⁶ PROX1 was also noted to be downregulated in familial combined hyperlipidemia, and it has been argued that this results in reduced clearance of plasma lipids in this disease.³⁷  It is possible that the identified variants in PROX1-AS1 increase the risk of osteonecrosis by altering lipid trafficking either within the compartment of the bone marrow or by increasing plasma lipids, a previously identified risk factor.¹⁰ 

The importance of the identified SNPs in osteonecrosis is further supported by the findings of the enhancer enrichment analysis. The top 92 validated SNPs were enriched for locations within enhancers active in tissues closely related to osteonecrosis, specifically mesenchymal progenitors (Table 2). These cells, under the influence of cytokines including BMP7³⁸  and environmental signals including glucocorticoids,³⁹  vitamin D, and mechanical load²⁹  can differentiate into either osteoblasts or adipocytes.^40,41  SNPs significantly associated with osteonecrosis (P < .001) were linked to 7 genes in the adipogenesis pathway (Table 3), supporting the importance of genes affecting mesenchymal differentiation in osteonecrosis. These findings are additionally supported by the findings in prior GWAS linking variants in fat and cholesterol metabolism genes (including ACP1) with osteonecrosis.¹⁰ 

Pathway analysis demonstrated that variation in glutamate receptor signaling contributes to osteonecrosis in children <10 years of age. Variants in 6 genes in the glutamate receptor pathway were associated with the development of osteonecrosis in this cohort, including the top validated nonsynonymous variant (rs34144324 in GRID2; Table 3). Interestingly, the glutamate receptor signaling pathway was the top pathway represented by genetic variants in a cohort of HR patients of all ages.¹⁵  Although glutamate receptor signaling appears to play a role in osteonecrosis across the age spectrum, the association appears stronger in older patients than in young patients, as both the top SNP and top pathway in a cohort dominated by older patients were related to glutamate signaling.¹⁵  In contrast, in young patients, genetic variation in genes expressed in mesenchymal and adipose tissues appear more important as supported by single SNP, enhancer, and pathway analyses (Tables 1-3). It is also possible that the apparent stronger association of glutamate receptor signaling in the prior study reflects the effects of high-risk therapy in addition to differences in age. However, prior findings in the Vanderbilt BioVU cohort, in which most patients were >10 years old and who developed osteonecrosis after glucocorticoids outside the treatment of ALL,¹⁵  suggests this association may be more due to age than therapy intensity.

Prior genome-wide studies^10,15  failed to identify an association between osteonecrosis and BMP7. This is likely due to the greater number of older patients in prior genome-wide studies. The reasons for the specificity of BMP7 to osteonecrosis risk in younger patients but not older patients are not clear. This may be explained by interactions between BMP7 and vitamin D, as vitamin D levels vary across the pediatric age range.^42,43 BMP7 expression in the absence of vitamin D induces osteoblast differentiation and mineralization, but this effect is reversed in the presence of 1-25(OH)₂ vitamin D₃.⁴⁴  Thus, it is possible that variation in BMP7 is important for osteonecrosis risk in younger but not older children due to age related differences in vitamin D levels. Alternatively, because BMP7 can induce angiogenesis,^45,46  it is possible variants in BMP7 have a smaller impact on osteonecrosis in older children and adolescents when bony growth is complete and the bony vasculature is already established.¹⁵  The BMP7 and PROX1-AS1 variants replicated in AALL0232 patients <10 years of age (P = .0498 and P = .008, respectively) but not in patients ≥10 years of age (P = .73 and P = .74, respectively), which also suggests the effects of the variants are age rather than treatment specific. Future evaluation of these findings in additional SR and high-risk patient populations under age 10 will allow further understanding of the interaction of younger age and treatment intensity in the development of osteonecrosis.

Although less common in children <10 years of age, osteonecrosis remains a major therapeutic toxicity in this age group. The findings of this GWAS provide new insight into the genetic features associated with the development of osteonecrosis in younger patients while verifying that glutamate receptor genetic variation may be important across all age groups.¹⁵  Our findings support the notion that some genetic risk factors for osteonecrosis are age specific, whereas others are not.

This study was supported by the National Institutes of Health, National Institute of General Medical Sciences (grants GM92666, GM115279, and GM97119), National Cancer Institute (grants CA21765, CA36401, CA142665, CA156449, CA98543 [COG Chair's grant], CA98413 [COG Statistical Center], and CA114766 [COG Specimen Banking]), Leukemia & Lymphoma Society (grant 6168), and the American Lebanese Syrian Associated Charities.

Contribution: E.A.R. and M.V.R. conceived and designed the study; S.E.K., L.A.M., W.Y., K.W.M., C.S., C.L., T.Y.C., M.L.L., S.P.H., M.D., and M.V.R. collected and assembled the data; S.E.K., W.Y., C.S., T.Y.C., C.C., and M.V.R. analyzed and interpreted the 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 Place, Room I-5112, Memphis, TN 38105; e-mail: mary.relling@stjude.org.