|Home | About | Journals | Submit | Contact Us | Français|
Treatment with glucocorticoids is associated with lower bone mineral density (BMD). We performed a genome-wide association study to analyze interactive effects between genotypes and cumulative dose of prednisone (PD) over 4.3 years of follow-up period on the final BMD Z-scores in 461 white children from the Childhood Asthma Management Program. No variants met the conventional criteria for genome-wide significance, and thus we looked for evidence of replication. The top 100-ranked single nucleotide polymorphisms (SNPs) were then carried forward replication in 59 children with acute lymphoblastic leukemia (ALL) exposed to large fixed doses of PD as part of their chemotherapeutic regimen. Among them, rs6461639 (interaction P = 1.88 ×10−5 in the CAMP population) showed a significant association with the final BMD Z-scores in the ALL population (P = 0.016). The association of the ALL population was only present after correction for the anti-metabolite treatment arm (high vs. low dose). We have identified a novel SNP, rs6461639, showing a significant effect on the final BMD Z-scores in two independent pediatric populations after long-term high dose PD treatment.
Glucocorticoids (GCs) are the core therapeutic agents in many childhood diseases. However, GCs have adverse effects on bone health as a consequence of decreased bone formation and increased bone resorption.1 It is known that long term use of GCs is associated with lower bone mineral density (BMD) compared to the population mean in children with nephrotic syndrome,2 systemic lupus erythematosus,3 acute lymphoblastic leukemia (ALL),4 and asthma.5 With ongoing skeletal development, children are more susceptible than adults to GCs effects on BMD. Various factors including gender, age, height, and pubertal stage are known to affect bone loss associated with GCs use in children.6 In terms of genetic factors, we previously reported that CRHR1 polymorphisms impacts the risk of bone density deficits in ALL children treated with GCs and antimetabolites.7 To date, there has been no study evaluating a genetic influence on GCs-induced bone loss in children using a genome-wide approach. We therefore performed a genome-wide association study (GWAS) with the hypothesis that we could identify novel genetic markers showing interactive effects with cumulative dose of GCs on the final BMD in children after a long-term GCs treatment. First, we tested our hypothesis by conducting a GWAS in white children with asthma participating in the Childhood Asthma Management Program (CAMP) trial.8 We then tested associations of the top 100-ranked SNPs in long-term survivors of childhood ALL who participated in our previous study.7
Each study was approved by the Institutional Review Board of the corresponding institution and informed consents were obtained from all study participants.
The primary group of subjects consisted of Caucasian children from the CAMP trial. The demographics of them and the study design have been described elsewhere.8 In brief, 1,041 children with mild-to-moderate asthma aged 5 to 12 years were randomized to budesonide, nedocromil, or placebo and followed for a mean of 4.3 years. For asthma exacerbations, short courses of oral prednisone (PD) were prescribed per protocol. Each burst consisted of 2 mg/kg per day (up to 60 mg) of PD for 2 days followed by 1 mg/kg per day (up to 30 mg) for 2 days and an option to continue dosing was allowed if improvement was insufficient. The CAMP population was divided into four groups according to their cumulative dose of PD over the follow-up period; PD non-takers and tertiles of PD takers [1st tertile (range: 5 ~ 330mg), 2nd tertile (range: 340 ~ 870mg) and 3rd tertile (range: 880 ~ 4,660mg)]. BMD measurements (in grams per square centimeter) of the lumbar spine (L1–L4) were performed yearly during the study period by means of dual-energy x-ray absorptiometry with the Hologic (Waltham, Mass) QDR-1500 (at 6 centers) or the Lunar (Madison, Wis) DPX (at 2 centers). Hologic DEXA machines were further divided by the use of pencilbeam or fan-beam measurements. Detailed methods to convert Lunar measures to Hologic values and to adjust deviations between pencil- and fan-beam measurements were described in our previous report.9 The BMD Z-score that represents the difference between the measured bone density in a subject and the average bone density in age- and sex-matched controls was calculated by using the CAMP internal references.5 In the present study the final BMD Z-score measured after mean follow-up duration of 4.3 years was used as a target phenotype. We confirmed that the chance of delayed puberty which might affect the final BMD Z-scores was small in the CMAP subjects.
The secondary population was drawn from participants in our previous study with ALL children.7 They were treated at St Jude Children’s Research Hospital between 1994 and 1998.10 For induction therapy, PD 40 mg per square meter of body surface area per day (mg/m2/day) was administered for one month, thereafter patients received dexamethasone 8 mg/m2/day (roughly 64 mg/m2/day in PD equivalents) for 7 days every 4 weeks for 120 weeks of continuation therapy. Based on the probability of ALL relapse, they were divided into two treatment arms; high-dose anti-metabolites + GC vs. low-dose anti-metabolites + GC. The high exposure to anti-metabolites could have obscured genetic effects on GCs-induced bone loss.7,11,12 We hypothesized that different levels of anti-metabolite administered might influence the association of GC with the final BMD, although GC dose was fixed. The final BMD measurements of the lumbar spine (L1-L2) were performed by Quantitative Computed Tomography with a Siemens Somatom-Plus spiral CT scanner (Siemens, Iselin, NY) at least 4 years after completing ALL therapy. Z-scores of the final BMD (expressed as the mean of the L1-L2 bone) were calculated for each patient. Detailed methods were described in our previous report.7
The CAMP population was genotyped on the HumanHap550v3 BeadChip or Infinium HD Human610-Quad BeadChip by Illumina, Inc. (San Diego, CA) and the ALL population was genotyped on the Affymetrix 5.0 chip by Affymetrix, Inc. (Santa Clara, CA). Marker quality control is detailed in the online supplement. Each single nucleotide polymorphism (SNP) had a completion rate > 95%, a minor allele frequency (MAF) of > 0.05, and a Hardy-Weinberg equilibrium P value of > 0.001. After merger with the phenotype dataset, 461 subjects with asthma and 59 subjects with ALL were included in the final analysis.
Using a linear regression gene-environment interaction model as implemented in the PLINK (http://pngu.mgh.harvard.edu/purcell/plink/),13 we performed a GWAS to test the following model: final BMD ~ genotype (additive genetic model) + cumulative dose of PD + genotype*cumulative dose of PD + covariates. Based on our previous report, age, sex, BMD (Z-score), height, body mass index, serum concentration of vitamin D, and Tanner stage measured at baseline were included as covariates.9 Since we focused on the interactive effect between genotype and cumulative dose of PD on final BMD, we reported P values from the ‘genotype*cumulative dose of PD’ term. We then selected top-ranked 100 SNPs to carry forward to test whether significant associations were found in the ALL population. Analysis was done using the following model: final BMD ~ genotype (additive genetic model) + gender + obesity + ancestry + treatment arm (high dose anti-metabolites group vs. low dose anti-metabolites group). A significant association in the secondary population was defined as a result with a nominal P value of less than 0.05. We confirmed that the assumptions of the linear regression model were met. Analysis was performed with R version 2.15.3 (www.r-project.org).
Table 1 summarizes the baseline characteristics of the primary (CAMP) and the secondary (ALL) populations. The final BMD Z-score of the third tertile of PD takers was significantly lower than that of PD non-takers (P = 0.015, Supplementary Figure 1).
No variants met the conventional criteria for genome-wide significance, and thus we looked for evidence of replication. The top 100-ranked SNPs identified in the primary GWAS are listed in Table S1 (Supplementary Table 1). Among those SNPs, we found that rs6461639 (P = 1.88 ×10−5) also showed a significant effect on the final BMD in the ALL population (P = 0.016). Power for the interaction analysis was calculated using Quanto software 1.2.4 (http://biostats.usc.edu/Quanto.html). Quanto is a program for computing either power or required sample size for association studies of genes, environmental factors, gene-environment interaction.14 An additive genetic model and alpha value (two-sided) of 0.05 were assumed. For rs6461639, an estimated power was 0.39 for main (gene) effect and 0.30 (2 df test) for interaction. The CAMP subjects carrying two variant (minor) alleles of rs6461639 showed significantly lower final BMD Z-scores compared to those with two reference (major) alleles after adjustment by covariates [P = 0.0016, odds ratio (95% confidence interval) = 0.54 (0.37 – 0.78), Figure 1a] in the third tertile of PD takers. The same trend was found in the ALL population (Figure 1b). Given that dose of PD administered in the ALL population was much higher than that of the CAMP population, we may say that the ALL population was in similar situation with the third tertile of PD takers (higher cumulative dose of PD) in the CAMP population.
In addition to tertiles analysis, we assessed the linear association (e.g. ’final BMD Z-score = a + b * cumulative dose of PD (mg) (a; intercept, b; slope)’) for each genotype of rs6461639 in the CAMP population (Figure 2). The significant linear correlation was seen only in subjects with two variant alleles of rs6461639 (Pearson’s correlation coefficient = −0.322, P = 0.0074). To further explore potential clinical effects behind our association, we assessed the relationship between non-PD takers and a hypothetical subject with asthma who was homozygous for the mutant rs6461639 allele following 4,000 mg of corticosteroid use over time (Figure 3). Overall, this subject would demonstrate −1.25 as the final BMD Z-score of the CAMP subject vs. the non-PD counterpart.
In the present study, we identified top-ranked 100 SNPs that showed significant interactive effects with cumulative dose of oral PD on the final BMD Z-score in childhood asthmatics. Among those 100 SNPs, we confirmed that rs6461639 located on the intron of Rap guanine nucleotide exchange factor 5 (RAPGEF5) gene was also significantly associated with the final BMD Z-score in ALL patients treated with long-term PD after adjustment of anti-metabolite effects. The significant interaction in the CAMP cohort was mainly driven by the effect observed in high dose PD takers (tertile 3 in Figure 1a). Considering that the children with ALL were treated with much higher doses of PD, the ALL cohort provides a possible evidence for replication of results shown in children taking high doses of oral GCs in the CAMP cohort. The pharmacogenomic magnitude of the clinical effect is substantial: the final BMD Z-score would be 1.25 less in CAMP subjects with two variant alleles of rs6461639, if we assumed that total dose of PD during 4.3 years of follow-up was 4,000 mg and factors other than genotype did not affect BMD. A total dose of 4,000 mg is equivalent to 5 mg of PD every other day for 4.3 years, which is not unusual in the real world setting. A Z-score of −1.25 indicates that a subject’s BMD would be within the lowest 10.6% for age. To the best of our knowledge, this is the first study showing an interactive effect on BMD between genetic factors and cumulative dose of GCs.
Bone mineralization throughout childhood and adolescence determines the peak bone mass. It is thought that heritable factors account for an estimated 60–80% of the variability in peak bone mass, with diet, physical activity and hormonal status serving as important modifiers of bone accrual.15 As discussed earlier, a chronic administration of GCs evidently hindered proper bone mineralization in children.2–6 GCs-induced bone loss is known to be associated with combined inhibition of both the osteoblast and osteoclast indicating that GCs-mediated inhibition of osteoclast function retards remodeling and thus dampens osteogenesis as shown in other states of the dynamic bone.16,17 The issue as to how GCs affect bone resorption by osteoclast has been controversial. A recent report comparing the effects of dexamethasone on wild-type osteoclasts with those derived from mice with disruption of the GC receptor in osteoclast lineage cells showed that the inhibitory effect of dexamethasone on bone resorption reflected failure of osteoclasts to organize their cytoskeleton.17 The osteoclast undergoes continuous reorganization with different phases of the resorptive cycle and thus a unique cytoskeleton is one of important characteristics of osteoclast.18 Dexamethasone in low nanomolar concentrations disarranged the osteoclast cytoskeleton, resulting in failure to spread and ineffectively resorb mineralized matrix.17 The osteoclast cytoskeleton is modulated by small GTPases like RhoA19 and Rac20 that transit to their GTP bond state under the influence of guanine nucleotide exchange factors (GEFs). Therefore, it is possible that GEFs act in the pathogenesis of GCs-induced decrease in BMD. Rs6461639 is located on the intron of RAPGEF5 gene that serves as Rap1 activators by promoting acquisition of GTP to maintain the active GTP-bound state.21 Notably, osteoclast-specific deletion of Rap1 yielded osteopetrotic mice by promoting talin/β integrin recognition and thus changing cytoskeletal organization.22 In addition, PPI finder, a web-based mining tool for human protein-protein interactions23 indicates that RAPGEF5 has a GC-regulated promoter.24 Taken together, RAPGEF5 regulated by GCs may deteriorate BMD by deranging cytoskeletal organization of osteoclast via Rap1. Recent genome-wide meta-analysis have identified several loci associated with BMD including 18p11.21 (FAM210A), 7q21.3 (SLC25A13), 11q13.2 (LRP5), 4q22.1 (MEPE), 2p16.2 (SPTBN1), 10q21.1(DKK1), 1p36.12 (ZBTB40), 1p31.3 (GPR177), and 7q.31.31 (CPED1 and WNT16).25–27 However, studies evaluating BMD decreases induced by GCs were not included, which may be the reason why RAPGEF5 was not listed in the previous results.
Rs6461639 is located outside of the coding region of RAPGEF5 gene, so we search its regulatory potential using ‘RegulomeDB’ which integrates a large collection of regulatory information from high-throughput experimental data sets and provides functional categories based on the strength of evidence.28 Rs64614639 is classified as the ‘likely to affect binding’ category (2b) based on direct evidence of binding through ChIP-seq and DNase with either a matched positional weight matrix to the ChIP-seq factor or a DNase footprint.29 We have identified 2 SNPs (rs1981597 and rs6968584) in LD with rs6461639 using the SNAP site (r2 threshold = 0.8, distant limit = 500 kb, CEU population).30 However we failed to find any functional relevancy of those two SNPs. In addition, we have identified 35 SNPs (on the same chromosome, frequency > 0.1, but not rs6461639 or in significant LD with rs6461639) showing significant (P < 0.0001, CEU population) effects on the expression of RAPGEF5 using the SCAN database.31 Although we cannot confirm a direct link between rs6461639 and RAPGEF5 gene function, rs6461639 is a functionally relevant SNP and thus may have its own role in the pathogenesis of GCs-induced decrease in BMD. Rs6461639 was significantly associated with the final BMD Z-score only in childhood asthmatics with a median cumulative dose of PD greater than 332.5 mg per year; this result was replicated in childhood ALL patients treated with cumulative dose of PD much higher than childhood asthmatics. It is not uncommon to see asthmatic children (and adults) who require equivalent doses of GC to control their disease activities in real practice. Considering that rs6461639 is a common SNP with MAF of 34% in Caucasians,32 our findings suggest that there may be clinical utility in genotyping rs6461639 in children before initiating GCs treatment, especially in cases of high dose therapy, to effectively anticipate and manage the possible decrease in BMD.
There are several potential limitations in generalizing our results. First, a P value of rs6461639 does not exceed the conventional levels for genome-wide significance (P < 5.0 × 10−8) in the first cohort and for Bonferroni correction adjusted P value (P = 0.005) in the replication cohort. Therefore, we cannot completely exclude the possibility that our results may be due to chance alone (false-positive associations). One possible explanation for these P values is the small sample size of the primary cohort was one of the contributing factors, as large sample sizes are generally needed when interactions are introduced into hypothesis testing.33 In addition, it is possible that the current GWAS design does not provide enough statistical power to detect interactions, because such a design mainly focuses on the main effect of a signal SNP on risk of common, complex diseases.34 However, it was noted that even extremely small P values provide surprisingly little information as the likelihood of subsequent replication.35 For examples, genetic variations with genome-wide significance showed poor consistency across GWAS of major depression36 and allergic rhinitis.37 Hence, we looked for evidence of replication of SNPs from our first GWAS. Secondly, the ALL population was not ideal to reveal interactive effects of SNP and the cumulative dose of GCs, because they were treated by only fixed dose of GCs, although we assumed that different levels of anti-metabolite might influence the association of GCs with the final BMD. However, it was very difficult for us to find clinical cohorts taking moderate to large doses of oral GCs with available BMD and GWAS genotype data. Therefore, it should be acknowledged that the association of the ALL population was only present after correction for the anti-metabolite treatment arm. Thirdly, the time gap between the measurement of clinical variables at baseline (covariates in our analysis) and the measurements of final BMD (target phenotype in our analysis) was considerable (over 4 years) and we could not measure the possible changes of all clinical variables that may have affected BMD during follow-up period. Nevertheless, we have adjusted our analysis for well-known potential confounders related to BMD, including age, sex, baseline BMD (Z-score), height, body mass index, serum concentration of vitamin D, and Tanner stage measured at baseline. Finally, although we found that rs64616395 has a regulatory potential, further mechanistic studies are needed to confirm the effect of rs6461639 on RAPGEF5 gene.
In conclusion, we have identified rs6461639 which demonstrates a pharmacogenomic interactive effect with cumulative dose of PD on the final BMD Z-score in children after a long-term GCs treatment. These results widen our understanding of the mechanistic basis as to how GCs induce decreases in BMD and may eventually be used as a preventive strategy in practice.
This work was supported by the NIH Pharmacogenomics Research Network: U01 HL65899 and U01 GM92666, and by NIH NHLBI R01 HL092197, NINR R01 NR013391, and NCI grants CA 142665, and CA 21765, and by the American Lebanese Syrian Associated Charities (ALSAC). The Childhood Asthma Management Program is supported by contracts NO1-HR-16044, 16045, 16046, 16047, 16048, 16049, 16050, 16051, and 16052 with the National Heart, Lung, and Blood Institute and General Clinical Research Center grants M01RR00051, M01RR0099718-24, M01RR02719-14, and RR00036 from the National Center for Research Resources.
CONFLICT OF INTEREST
All the authors state that there are no conflicts of interest to disclose.
The authors directed and had access to all the analyses and the full clinical and genetic database, wrote all drafts of the report, decided to publish the results and attest for the accuracy and completeness of the data.
Supplementary information is available at The Pharmacogenomics Journal’s website