A coding variant in RARG confers susceptibility to anthracycline-induced cardiotoxicity in childhood cancer.

Anthracyclines are used in over 50% of childhood cancer treatment protocols, but their clinical usefulness is limited by anthracycline-induced cardiotoxicity (ACT) manifesting as asymptomatic cardiac dysfunction and congestive heart failure in up to 57% and 16% of patients, respectively. Candidate gene studies have reported genetic associations with ACT, but these studies have in general lacked robust patient numbers, independent replication or functional validation. Thus, the individual variability in ACT susceptibility remains largely unexplained. We performed a genome-wide association study in 280 patients of European ancestry treated for childhood cancer, with independent replication in similarly treated cohorts of 96 European and 80 non-European patients. We identified a nonsynonymous variant (rs2229774, p.Ser427Leu) in RARG highly associated with ACT (P = 5.9 × 10(-8), odds ratio (95% confidence interval) = 4.7 (2.7-8.3)). This variant alters RARG function, leading to derepression of the key ACT genetic determinant Top2b, and provides new insight into the pathophysiology of this severe adverse drug reaction.

Introductory

Anthracyclines are used in over 50% of childhood cancer treatments[1], but their clinical utility is limited by anthracycline-induced cardiotoxicity (ACT) manifesting as asymptomatic cardiac dysfunction and congestive heart failure in up to 57% and 16% of patients, respectively[2,3]. Candidate gene studies have reported genetic associations with ACT[4-22] but in general they lack robust patient numbers, independent replication or functional validation. Thus, individual ACT susceptibility remains largely unexplained[12,13]. We performed a genome-wide association study in 280 patients of European ancestry treated for childhood cancer with independent replication in similarly-treated cohorts of 96 European and 80 non-European patients. A nonsynonymous variant (rs2229774, p.Ser427Leu) in RARG was highly associated with ACT (P=5.9×10−8, odds ratio (OR) (95% confidence interval) = 4.7 (2.7–8.3)). This variant alters RARG function leading to derepression of the key ACT genetic determinant, Top2b, and provides novel insight into the pathophysiology of this severe adverse drug reaction (ADR).

The etiology and pathogenesis of ACT are poorly understood. Several clinical[23-25] and cytotoxicity[26,27] studies support a strong genetic contribution. The identification of genetic risk factors for ACT contributes to a better understanding of its pathophysiology, which could lead to new approaches to predict, prevent and treat this severe ADR[12,13].

Accordingly, we recruited well-phenotyped patients treated with anthracyclines for childhood cancer (≤18 years at treatment) from 13 pediatric oncology centers across Canada. Clinical information was used to assess cardiac function, define cases and controls, and determine important baseline differences between these groups. Cases were defined as exhibiting shortening fractions [SF] of ≤24% or signs and symptoms of cardiac compromise requiring intervention based on CTCAEv3, while controls had SF ≥30% and no symptoms of cardiac compromise for at least 5 years after treatment[12,13]. Detailed methodological descriptions are available in the Online Methods.

A stage 1 discovery analysis (Supplementary Figure 1) was performed in Canadian patients of European ancestry (280 patients; 32 cases and 248 controls). Compared with controls, cases were significantly older at the start of treatment, had higher cumulative anthracycline exposure (dose), were less likely to have acute lymphoblastic leukemia (ALL), but more likely to have rhabdomyosarcoma or Ewing’s sarcoma, and receive radiotherapy to the heart (RT) (P<0.05; Supplementary Table 1). Age, dose and RT are established risk factors for ACT[28,29] and after verifying that the associations with ALL, Ewing’s sarcoma and rhabdomyosarcoma tumor types were not exclusively explained by the cumulative anthracycline dose (remained associated with ACT after accounting for dose), these six clinical factors were included as covariates for logistic regression analyses. Notably, the limited sample size may have precluded our detection of other significant clinical differences between cases and controls.

Patients were genotyped using an Illumina HumanOmniExpress (740K SNP) assay with 657,694 SNPs passing quality control assessment, conferring a Bonferroni-adjusted multiple testing threshold of P<7.6×10−8 (Online Methods). The GWAS discovery analysis was performed using logistic regression adjusted for age, dose, ALL, Ewing’s sarcoma, rhabdomyosarcoma, and RT (Supplementary Figure 2a). Analysis of the test statistics (λGC = 1.021), suggested that they were not influenced by cryptic population stratification (Supplementary Figure 2b). Eighteen variants tagging nine distinct linkage disequilibrium (LD) blocks (r2>0.9 and D′>0.9 in the 1000 Genomes CEU reference dataset) with P<1.0×10−5 were prioritized for further analysis in an independent patient population (Supplementary Table 2).

In stage 2, we genotyped one GWAS candidate variant per LD block in an independent Dutch population of childhood cancer patients of European ancestry (96 patients; 22 cases and 74 controls). In this cohort cases and controls were similarly matched with the exception that cases had higher cumulative anthracycline exposure (Supplementary Table 1). Genetic associations were tested using logistic regression with adjustment for cumulative anthracycline dose. Of the 9 candidate variants tested, only rs2229774, a nonsynonymous coding variant (p.Ser427Leu) in RARG (Retinoic Acid Receptor Gamma), was replicated (P=0.0043, OR=4.1 (1.5–11.5); Fig. 1, Table 1).

Figure 1

A pharmacogenetic association with susceptibility to anthracycline-induced cardiotoxicity is situated within RARG

a) Association results are shown for genotyped (circles) and imputed (squares) SNPs along with recombination rates for a 122 kb region of chromosome 12q13.13. Each point represents the nominal P-value (left y-axis) for the stage 1 cohort. P-values are from logistic regression analysis using an additive model, adjusted for age, cumulative dose, tumor type (ALL, Ewing’s sarcoma and rhabdomyosarcoma) and cardiac radiation therapy. SNPs are colored according to their pairwise correlation (r2) with rs2229774 (purple circle) using the 1000 Genomes CEU reference population. Overlaid are the recombination rates (right y-axis) for estimating putative recombination hotspots also based upon the 1000 Genomes CEU population. b) The linkage disequilibrium (D′) based upon the 1000 Genomes CEU population depicted for this region similarly demonstrates the associated haplotype is localized to RARG. Details of D′ color coding are provided in Online Methods.

Table 1

Association of RARG rs2229774 with Anthracycline-induced Cardiotoxicity in Childhood Cancer Patients

Biomarker			Pharmacogenomic Analyses			Adjusted Logistic Regression (Additive Model)		Genotypic Test
SNP	Gene	Function	Study Population	MAF Cases	MAF Controls	P	Odds Ratio (95%CI)	P
rs2229774	RARG	NON-SYN CODING	Stage 1 – Discovery GWASa Canadian European Patient	0.297	0.081	5.0 × 10−6	7.0 (2.9 – 17.0)	4.1 × 10−8
			Stage 2 – Replicationb Dutch European Patients	0.25	0.061	0.0043	4.1 (1.5 – 11.5)	0.0042
			All European Patientsa	0.278	0.076	7.8 × 10−8	5.4 (2.9 – 10.3)	1.2 × 10−9
			Stage 3 – Replication Non-European Patients	0.158	0	N/Ac	N/A	1.2×10−4
			All Populationsa European and Non-European Patients	0.247	0.064	5.9 × 10−8	4.7 (2.7 – 8.3)	4.3×10−11

Covariates for the Logistic regression were age at treatment, cumulative anthracycline exposure, radiotherapy involving the heart and incidence of acute lymphoblastic leukemia, rhabdomyosarcoma and Ewing’s sarcoma.

Covariate for the Logistic regression was cumulative anthracycline exposure.

Not applicable, rs2229774 absent in controls.

Given that the Stage 1 and 2 cohorts differed in some patient characteristics, such as anthracycline dose, the replication analysis merits circumspect interpretation. To address this we performed additional genetic association analyses by logistic regression with the prioritized GWAS variants in subsets of cases and controls stratified by low-to-moderate (≤250 mg/m2) and high (>250 mg/m2) anthracycline exposure (Supplementary Table 3). RARG rs2229774 was associated with ACT in stages 1, 2 and the combined analysis at both low-to-moderate (P=4.1×10−4, P=0.0036, and P=9.8×10−6, respectively) and high anthracycline doses (P=0.0021, P=0.084, and P=8.7×10−4, respectively). By contrast, none of the other top GWAS candidate SNPs were significantly associated with ACT in the replication cohort, even at low-to-moderate anthracycline exposure. Logistic regression analyses in the Stage 1 cohort indicated that RARG rs2229774 was similarly associated with early-onset chronic ACT (n=16 cases; P=2.8×10−4; OR=7.3 (2.3–22.9)) and late-onset chronic ACT (n=16 cases; P=5.0×10−3; OR=5.8 (1.7–19.4)).

In stage 3, we examined the association of RARG rs2229774 with ACT in a third cohort of non-European patients representative of different ancestries (80 patients; 19 cases and 61 controls; Supplementary Table 4). Due to the absence of RARG rs2229774 in controls, logistic regression with adjustment for population stratification and other clinical covariates could not be performed in this cohort (Supplementary Table 5). Instead the association of RARG rs2229774 with ACT in this cohort was analyzed by genotypic test. RARG rs2229774 was highly associated with ACT (P=1.2×10−4; Table 1) in the stage 3 cohort, and in each of the four non-European populations separately (African, Aboriginal Canadians, Hispanic, and East Asian; Supplementary Table 5). Notably, genotypic association testing reached genome-wide significance in the discovery and combined cohorts (P<5.0×10−8; Table 1).

The stage 3 cohort, comprised of four distinct populations, was assessed for population stratification using the genomic inflation factor and principal component analysis (PCA). PCs 1–10 did not significantly differ between cases and controls (P>0.14) and the λGC of 0.941 suggested no confounding effect of population stratification on the genetic association of rs2229774 with ACT in this cohort. In line with the reported high incidence of ACT in studies conducted in India[30] and Pakistan[31], and in African-American patients in the USA[23,24], RARG rs2229774 is more frequent in South Asian (22%) and African populations (11%) compared to European (6%) and Hispanic (5%) populations. By contrast, rs2229774 is very rare in East Asian (0%) populations (1000 Genomes). It will be interesting to correlate the frequency of rs2229774 with the incidence of ACT in children of different ancestries as data on multiethnic studies become available for comparison.

A meta-analysis of all study populations (456 patients; 73 cases and 383 controls) using logistic regression adjusted for age, dose, ALL, RT, Ewing’s sarcoma and rhabdomyosarcoma showed that RARG rs2229774 was significantly associated with ACT (P=5.9×10−8, OR=4.7 (2.7–8.3)), which surpassed the Bonferroni-adjusted multiple testing correction threshold, with no significant heterogeneity across the GWAS and two independent replications (Pheterogeneity=0.402). Overall, rs2229774-carriers had significantly increased odds of developing ACT compared to non-carriers (OR=5.2 (3.0–9.0); P=5.9×10−10).

To fine-map the rs2229774 association with ACT, we imputed variants on chromosome 12 using the stage 1 cohort and the 1000 Genomes CEU reference population, and tested their association with ACT by logistic regression (Fig. 1). Of the 1,005,286 imputed SNPs, three were associated with ACT at P<1.0×10−5 (rs11170481, rs73309171 and rs57789211; Table 2). These were all located in introns of RARG, highly correlated with each other (D′≥0.91, r2≥0.54), and in high linkage disequilibrium (LD) with rs2229774 (D′=1.00, r2≥0.55). Since the ACT associations of the imputed variants were similar to rs2229774, and logistic regression with conditioning on rs2229774 abolished these associations (Table 2), the non-synonymous coding variant rs2229774 was prioritized for functional characterization to gain biological insight into its association with ACT.

Table 2

Fine Mapping of Genetic Association Signals in RARG gene regiona.

Biomarker				Pharmacogenomic Analyses			Adjusted Logistic regression (additive model)b
SNP rs-IDc	Positiond	Type	Source	LD with rs2229774 D′ (r2)e	MAF Cases	MAF Controls	P	Odds Ratio (95%CI)	Conditional Analysis on rs2229774 P
rs11170481	53611791	Intronic	Imputed	1.00 (0.84)	0.313	0.081	1.7×10−6	7.0 (3.0 – 16.6)	0.10
rs73309171	53606565	Intronic	Imputed	1.00 (0.55)	0.281	0.071	4.1×10−6	7.5 (3.0 – 18.4)	0.42
rs57789211	53609992	Intronic	Imputed	1.00 (0.84)	0.281	0.071	4.1×10−6	7.5 (3.0 – 18.4)	0.42
rs2229774	53605545	Nonsyn	Genotyped	1.00 (1.00)	0.297	0.081	5.0×10−6	7.0 (2.9 – 17.0)	–

1,005,286 additional variants on Chr12 imputed into stage 1 cohort using the CEU component of the 1000 Genomes populations as a reference.

Covariates for logistic regression analysis (additive model) were age at treatment, cumulative anthracycline exposure, radiotherapy involving the heart and incidence of acute lymphoblastic leukemia, rhabdomyosarcoma, Ewing’s sarcoma and rs2229774 where indicated.

Association analyses for imputed SNPs were restricted to those with BEAGLE allelic R2 ≥ 0.5.

Chromosomal positions in the GRCH37.p13.

Calculated using the CEU component of 1000 Genomes reference population.

Retinoic acid receptors bind to DNA regulatory sequences termed retinoic acid response elements (RARE) and transcriptionally co-regulate downstream gene expression in response to their agonist all-trans retinoic acid (ATRA)[32]. Notably, RARG can both activate and repress transcription in response to ATRA[33]. To explore the functional properties of rs2229774 (i.e. RARGS427L), we first examined its transactivation of RARG regulatory elements and then identified a putative role in the dysregulation of a critical gene involved in the development of ACT (see Online Methods). In HEK293T cells expressing RARGS427L, the ATRA-inducible transcriptional activation of a RARE-coupled luciferase reporter was significantly reduced compared to wild type RARG-expressing cells (Fig. 2a). Immunoblot analysis verified that both wild type and variant RARG proteins were detected at similar levels in HEK293T cell lysate (Fig. 2a inset). The significant 17% decrease in RARG activity conferred by the rs2229774 variant may be an underestimate since endogenous retinoic acid receptors were present in this assay. These results suggested that dysregulation of an RARG-regulated gene might underlie the association of rs2229774 with ACT.

Figure 2

Functional characterization of RARGS427L reveals impaired transcriptional regulation

a) Transcriptional activation of luciferase coupled to an optimized retinoic acid response element (RARE) by transiently transfected RARG wild type (WT) or RARGS427L in HEK293T cells. Averages of RARE-activation from aggregate data (n = 48; three independent experiments of sixteen replicates) are presented. Inset, Immunoblot of 20μg HEK293T lysate generated 48 hours post-transfection with empty vector (negative), or the indicated construct using anti-DDK 4C5 (top panel) and anti-GAPDH (bottom panel) antibodies. Untagged wild type RARG has an estimated molecular weight of 50.4 kDa. Molecular sizes are indicated on the left. b) Relative Top2b expression in untransfected or RARG-transfected H9c2 cells in the presence or absence of ATRA. Average values of relative Top2b expression from aggregate data (n = 12; four independent experiments of three replicates) are presented. (c) Repression of Top2b expression in RARG WT- or RARGS427L-transfected H9c2 cells compared to untransfected cells. Averages of normalized Top2b repression from aggregate data (n = 12; two independent experiments of six replicates) are presented. ** and *** denote P < 0.005 and P < 0.0001, respectively, using t-test (a and c) or one-way ANOVA with Tukey post-test analyses (b). Error bars; s.e.m.

Anthracyclines are mechanistically and genetically linked with Topoisomerase II beta (Top2b). Anthracyclines exert their anticancer activity by binding and inhibiting Topoisomerase II[34]. In addition, Top2b is necessary for the development of ACT in a murine model[35], while in a rat cardiomyoblast (H9c2) cell line, Top2b levels are decreased by the ACT cardioprotectant, dexrazoxane[36]. RARG expression has been reported as “particularly high” in the heart (Nuclear Receptor Signaling Atlas) and is highly induced in murine cardiac cells following cardiac injury[37]. Since RARG has been shown to bind to the Top2b promoter[38] it was one potential candidate of RARGS427L dysregulation. We found that Top2b expression in H9c2 cells was significantly decreased when human RARG was added, and this effect was further exacerbated with the addition of ATRA (Fig. 2b). By contrast, the RARGS427L variant did not repress Top2b expression as effectively as wild type RARG (Fig. 2c). Taken together these results demonstrate that RARG represses Top2b expression and that rat cardiomyoblasts carrying rs2229774 express higher basal levels of Top2b, consistent with an increased susceptibility to ACT.

We have identified a novel genetic biomarker for ACT, RARG rs2229774, using a three-stage genetic association study combined with biological functional analyses. This discovery, despite originating from a relatively small number of patients, is likely due to its large effect size. In general, sample size limitation remains a major challenge in pediatric cancer pharmacogenomic studies[39,40], where the need for a homogenous study population and well defined clinical phenotypes results in reduced numbers of study patients, particular the number of affected individuals. However, clinically relevant genetic markers of severe ADRs such as ACT, abacavir-induced hypersensitivity reaction, and carbamazepine-induced Stevens-Johnson syndrome, are expected to have large effect sizes and therefore, the number of patients, particularly the number of affected individuals, required to uncover these genetic associations from large-scale genome screens, has been shown to be in the range of 10–100 patients[41] – similar to the number in this study.

Our study was ≥ 80% powered to detect genome-wide associations of P≤5×10−8 with per-allele OR≥4.5 and MAF≥0.10. There are likely additional genetic contributions to ACT including previously reported associations[4-22] that were not uncovered in the current GWAS, potentially owing to the strict ACT case definition and resulting small patient numbers used in this study. For example, ABCC2/MRP2, HNMT and POR, although not reaching genome-wide significance in this study, exhibited strong trends of association (P<0.05; Supplementary Table 6) in the same direction of effect as previously reported[4,8,12,16,22]. This underscores the complementary nature of GWAS and candidate gene association studies in identifying variants associated with pharmacogenomic traits.

In agreement with in silico predictions, our in vitro studies demonstrate that rs2229774 causes a relatively tolerated amino acid substitution that results in a moderate, but significant, reduction in RARE transcriptional activation. We established a genetic interaction between RARG and Top2b, where expression of the latter is significantly repressed by RARG in the presence of ATRA. Notably, the treatment of acute promyelocytic leukemia with high cumulative doses of anthracycline and ATRA has resulted in a significantly lower incidence of ACT[42,43]. We further showed that rs2229774 causes derepression of Top2b in cardiomyoblasts, directly linking this variant with ACT. These data are consistent with a model where rs2229774 carriers have higher basal levels of TOP2B in cardiomyocytes, conferring increased susceptibility to cardiotoxicity when treated with anthracyclines. Nevertheless, despite rs2229774 encoding a non-synonymous variant, it is plausible that the rs2229774 haplotype confers increased susceptibility to ACT through a regulatory mechanism. To our knowledge a comprehensive analysis of RARG-regulated gene expression in cardiomyocytes has not been reported and it is possible that additional genes in cardiomyocytes may contribute to the development of ACT when dysregulated in RARG rs2229774 carriers. For example, ATRA-regulated gene expression in cardiomyocytes is important for cardiac development and the progression of cardiomyocyte hypertrophy[44-46], though the specific contribution of RARG to these processes is unknown.

The identification of RARG rs2229774, genetically and functionally linked to ACT, provides a clinical tool that may be used to predict genetic risk and improve ACT risk stratification. ACT is a clinically significant ADR, and carriers of rs2229774 have 5-fold increased odds of ACT. This novel finding merits further exploration of the role of RARG in ACT and of the clinical utility of RARG rs2229774 predictive testing to inform ACT risk assessment.

ONLINE METHODS

Study Design

Supplementary Figure 1 briefly outlines the study design. In stage 1, a total of 280 European Canadian patients (32 cases and 248 controls) were used as the discovery cohort and genotyped on the Illumina Infinium HumanOmniExpress panel (740K), to perform a GWAS for ACT (stage 1). Markers that reached P < 1.0×10−5 (ref. 47) in the discovery cohort were tested for replication in 96 European Dutch patients (22 cases and 74 controls) in stage 2. RARG rs2229774 was further tested for association with ACT in non-European populations in stage 3. We also performed a combined analysis for all European patients (stages 1 and 2) and an overall association test for rs2229774 using all populations (stages 1–3). In addition, we performed genetic fine mapping analyses of the associated region and functionally characterized rs2229774 in vitro.

Genetic ancestry was self-reported and ascertained by PCA using the EIGENSTRAT method, with the GWAS Illumina 740K (Canadian patients), Illumina 4.5K (Dutch patients) and Illumina 8K (USA patients) SNP genotype data sets.

Analysis of Clinical Data

This study was approved by the individual ethics committees/institutional review boards of the universities and institutions where patients were enrolled. Written informed consent/assent was obtained from patients/parents/legal guardians in accordance with the Helsinki Declaration as revised in 2008.

Patients’ medical records were reviewed prior to genotyping by a clinical pharmacologist, a cardiologist, an oncologist and an ADR surveillance clinician, who reviewed the echocardiogram test results and other clinical information. The demographic, clinical and therapeutic information was extracted from these medical records and included the following data: demographics, disease characteristics, chemotherapy, diagnostic echocardiograms to document baseline and follow-up cardiac function and any cardiac compromise and its severity and any symptoms and/or signs consistent with ACT. All patients were children ≤ 18 years at cancer diagnosis, who received anthracyclines as part of their chemotherapy protocol, had normal cardiac function before anthracycline chemotherapy and were recruited from outpatient clinics and inpatient units. ACT was monitored by echocardiograms according to the “The Children’s Oncology Group Long-Term Follow-Up Guidelines for Survivors of Childhood, Adolescent, and Young Adult Cancers”. The recommended frequency of echocardiograms is every year, every two years or every 5 years post-treatment depending on the age at treatment, radiation with potential impact to the heart and the cumulative anthracycline dose. Grading of ACT based on this detailed clinical characterization was performed using the Cancer Therapy Evaluation Program, Common Terminology Criteria for Adverse Events version 3 (CTCAEv3) as previously described[12,13]. Patients with serious ACT were defined as those with grade 2 or higher CTCAE impairment of cardiac function after treatment with anthracycline.

ACT was defined as early- or late-onset left ventricular dysfunction assessed by echocardiogram measurements using shortening fraction [SF]) and/or symptoms (dyspnea, orthopnea, and/or fatigue) and/or signs (edema, hepatomegaly, and/or rales) of cardiac compromise requiring intervention based on the CTCAEv3. Due to variability in echocardiographic measurements, a conservative threshold of SF ≤ 24% was used to define asymptomatic cases who developed ACT. All cases had grade 2 or higher ACT. Also, grade 2–4 ACT is the point at which anthracycline chemotherapy protocols recommend clinical intervention such as heart failure treatment, halting or reducing anthracycline doses and switching to alternate treatments. Early-onset chronic ACT was defined as developing less than 1 year after start of treatment[48,49] while late-onset chronic ACT was defined as developing more than 1 year after start of treatment[50,51] since ventricular dysfunction, heart failure, and arrhythmias can occur years or even decades after the discontinuation of anthracycline therapy[50,51]. To exclude transient acute ACT, only echocardiograms obtained ≥ 21 days after an anthracycline dose were considered. To exclude cardiotoxicity unrelated to anthracycline chemotherapy patients with no baseline echocardiogram were excluded from the study. Controls had no signs or symptoms of cardiac compromise at study participation (grade 0). Due to the delayed onset of ACT in some patients, SF ≥ 30% with ≥ 5 years of follow-up after the end of anthracycline treatment was used to define controls. The cumulative anthracycline exposure was calculated using the doxorubicin isotoxic equivalent[52] and cumulative dose stratification into low-to-moderate (≤ 250 mg/m2) and at high (> 250 mg/m2) anthracycline exposure was performed as previously described[6,14,53] where appropriate. Radiation therapy included significant radiation exposure to the heart or surrounding tissue. This included mantle and mediastinal radiation, whole lung radiation, whole or upper abdominal radiation, left sided flank radiation and total body irradiation.

Study Populations

Canadian Patient Populations

We genotyped 434 pediatric oncology patients treated with anthracyclines that were recruited from 13 pediatric oncology units from across Canada between February 2005 and April 2011 (refs. 12,13). A total of 13 samples failed quality control (details of quality control procedure described below). In 421 remaining patients (average call rate = 99.5%), we performed PCA with GWAS Illumina 740K SNP genotype data set to determine the population structure and defined four distinct population clusters comprised of European, African, East Asian and Aboriginal Canadian patients. Patients with genetic ancestry that fell outside of these four population clusters were excluded from further analysis. We excluded an additional 27 patients with CTCAE grade 1 toxicity (shortening fraction: 24% < SF < 30%). A total of 97 patients were excluded and 337 patients were available for further analyses: Europeans (280 patients; 32 cases and 248 controls), Africans (n = 11 patients; 2 cases and 9 controls), East Asians (n = 31 patients; 8 cases and 23 controls) and Aboriginal Canadians (n = 15 patients; 4 cases and 11 controls).

Dutch Patient Population

We recruited 128 pediatric oncology patients treated with anthracyclines from Emma Children’s Hospital/Academic Medical Centre in Amsterdam, the Netherlands, between July 2009 and April 2011 (refs. 12,13). We performed PCA using the Illumina 4.5K SNP genotype data set generated for these patients[54] and excluded 14 patients not genetically matching European ancestry. Based on the detailed clinical assessment of the 114 patients, we excluded another 18 patients with CTCAE grade 1 toxicity (shortening fraction: 24% < SF < 30%). A total of 96 patients (22 cases and 74 controls) were then available for further analysis.

USA Patient Population

We recruited 164 pediatric oncology patients treated with anthracyclines from Lucile Packard Children’s Hospital at Stanford (Palo Alto, USA) between December 2008 and October 2010. We performed a detailed clinical assessment of these patients and excluded 10 patients who had missing or abnormal baseline echocardiogram readings, 8 patients with CTCAE grade 1 toxicity (shortening fraction: 24% < SF < 30%), and 102 patients with insufficient cardiac follow-up data (either SF>30% but less than 5 year follow-up (97 patients), or SF<24% but <21 days from end of treatment (5 patients)). PCA identified four clusters of the remaining patients aligning with their self-reported ancestries of Hispanic, Asian, African-American, and European. The cluster of 23 Hispanic patients (5 cases and 18 controls) was selected for analysis in stage 3.

DNA Extraction and Molecular Genotyping

Genotyping experiments were conducted at the Canadian Pharmacogenomics Network for Drug Safety (CPNDS) genotyping core facility, Child & Family Research Institute, The University of British Columbia, Vancouver, BC, Canada. All patients provided a biologic specimen for DNA extraction. Genomic DNA was extracted using the QIAamp DNA purification system (Qiagen) and quantified by Quant-iT PicoGreen assay (Invitrogen), according to the manufacturer’s protocols. Whole Genome Amplification (WGA) was performed prior to genotyping for all samples available at a low concentration (< 50 ng/μl) and/or low volume (< 20 μl). The laboratory assistants were blinded to the case-control status of the patients genotyped in the study. To ensure the accuracy of all genotyping results, multiple positive and negative controls and replicate samples were included in all genotyping assays and plates. The concordance of genotype calls between replicate genotyped samples was 100%.

The genome-wide association study (GWAS) was performed using a customized Illumina Infinium HumanOmniExpress assay containing 738,432 SNPs according to the manufacturer’s instructions (Illumina). Genotypes were called with the Illumina Genome Studio software package and the SNPs were clustered using the Illumina 740K cluster file. A detailed GWAS quality control procedure was performed for all SNPs and samples prior to analysis using the Illumina Genome Studio software package.

The Dutch and Hispanic patients were genotyped for the GWAS candidate variants using TaqMan SNP genotyping assays (Applied Biosystems), according to the manufacturer’s protocols. The top replicated SNP from the Illumina Infinium HumanOmniExpress assay was re-genotyped in 100 randomly selected patients from the discovery patient population using the Taqman SNP genotyping assay. The concordance rate of genotype calls between the two assays was 100%.

GWAS quality control (QC)

The raw data was initially imported and clustered using the Illumina 740K cluster file. Next, an iterative genotyping/cluster QC process was performed using a sequence of QC filters for both SNPs and samples. Genetic markers were evaluated using a combination of thresholds for various quality control metrics available in GenomeStudio together with visual inspection of cluster plots for markers at the boundaries of the thresholds. Markers with call rates ≥ 95% were retained and poorly-clustered markers (call rates < 95%) were filtered out and re-clustered and the newly defined cluster positions were either left alone, manually edited, or dropped altogether after evaluating the call rates, cluster separation, mean normalized intensity, proximity of heterozygote clusters to a homozygote cluster, heterozygous excess, false homozygozity and reproducibility/replication errors. Deviation of the genotype distributions from Hardy-Weinberg equilibrium (HWE) was tested in control patients. All SNPs with Fisher’s exact test for HWE P-value < 1.0×10−4 were excluded.

We implemented the same quality-control procedures such as call rates and HWE from the GWAS, in the replication cohorts. Nine SNPs in the replication studies had call rates of >90%, while one SNP had a call rate < 90% and was subsequently excluded from further analysis.

Quality control for DNA samples was performed with SVS/HelixTree 8.1.1. Samples were excluded if they had a call rate below 95% and if the ancestry departed from the expected homogenous genetic ancestry. A total of 657,694 SNPs from the GWAS were retained after QC. Cluster plots for all GWAS associated SNPs were visually inspected.

Functional Assays

Constructs, cells and reagents

A Myc-DDK-tagged mammalian expression clone of human RARG was purchased from Origene. The rs2229774 SNP was introduced into this expression vector using the Quikchange II kit according to the manufacturer’s specifications (Agilent Technologies) with mutagenic primers (Supplementary Table 7). Transfections were performed with X-tremeGENE 9 (Roche) or Effectene (Qiagen) reagents according to the manufacturer’s specifications. HEK293T and H9c2(2-1) (ATCC CRL-1446) cells were purchased from ATCC (Cedarlane) and routinely cultured in DMEM supplemented with 10% FBS, 100U/ml Penicillin, 100μg/ml Streptomycin, with additional 0.25μg/ml Amphotericin B for H9c2 cells. RARG transcriptional regulation was assayed using the Cignal RARE Reporter luciferase kit (Qiagen). Luminescence assays were developed using the Dual-Glo Luciferase Assay System according to the manufacturer’s specifications (Promega) and read on a POLARstar Omega plate reader (BMG Labtech). For real-time RT-PCR (qPCR) experiments RNA was purified from H9c2 cells using the Ambion Purelink RNA mini kit with Purelink homogenizers and Purelink on-column DNAse digestion according to the manufacturer’s specifications (Life Technologies). cDNA was generated using the Invitrogen Superscript III first strand synthesis kit according to the manufacturer’s specifications (Life Technologies). Rat Top2b, rat Hprt1 and human RARG gene expression was measured using validated Applied Biosystems TaqMan assays Rn01537914_m1, Rn01527840_ml, and Hs01559234_m1, respectively (Life Technologies). qPCR was performed on the PikoReal 96 Real-Time PCR system (Thermo Scientific) and relative gene expression was calculated by the ΔΔCq method using the instrument software. The anti-DDK 4C5 antibody (catalog # TA50011) was purchased from Origene and the anti-GAPDH antibody (catalog # MAB374) was purchased from Millipore. Antibodies were used at (1:1,000) dilution.

RarG transcriptional regulation assay

To assess whether general RARG activity was affected by the rs2229774 variant, we used a reporter construct fused to an optimized RARE element (Cignal Reporter System, SABiosciences). HEK293T cells were reverse co-transfected in a 96-well plate with 50ng (per well) of either empty vector (pcDNA3.1), RARG WT or RARGS427L expression constructs and the Cignal RARE reporter Luciferase kit constructs according to the manufacturer’s specifications. After 20–24 hours of transfection cells were washed with PBS and the medium was replaced with 75μl Cignal Assay medium (OptiMem supplemented with 1 % charcoal-stripped FBS, 0.1mM NEAA, 1mM sodium pyruvate, 100U/ml Penicillin and 100μg/ml Streptomycin) containing 1μM ATRA or DMSO control for 6 hours. Firefly and Renilla luminescence was measured as indicated above. Firefly luciferase/Renilla luciferase ratios (L/R) were calculated for each well and converted to a Relative response ratio (RRR) = [(L/R)sample − (L/R)negative control]/[(L/R)positive control − (L/R)negative control] to allow sample comparison between experiments. The fold induction of ATRA-treated versus untreated samples was calculated using the corrected RRR values (subtraction of empty vector RRR values) for RARG WT- and RARGS427L-transfected samples.

Relative gene expression studies in rat cardiomyoblasts

7.5 × 104 H9c2 cells were seeded into each well of a 6-well dish in DMEM supplemented with 10% FBS. The following day, cells were transfected with 1μg RARG WT or RARGS427L expression constructs for 7 hours before replacing with fresh medium. Where required 250nM ATRA was added to cells 24 hours post-transfection. Cells were grown for 48 hours post-transfection then total RNA was immediately purified and used for cDNA synthesis. qPCR reactions were conducted in a 10μl reaction volume with standard cycling conditions according to the manufacturer’s specifications. Relative gene expression was calculated using Hprt1 as a housekeeping gene. To calculate the fold repression of Top2b expression in RARG-transfected H9c2 cells, relative Top2b expression in untransfected cells was divided by expression levels in transfected cells and normalized to the relative expression of the appropriate RARG construct.

Statistical Methods

We performed statistical analyses using SVS/HelixTree 8.1.1 (Golden Helix), R 3.1.0 for Statistical Computing, SPSS Version 18.0 (IBM), Quanto, Haploview[55], LocusZoom[56], Epi Info™ Version 7.1.3, Comprehensive Meta-Analysis, BEAGLE 4 (ref. 57), GraphPad Prism 5.0a and PikoReal version 2.1 software packages. All statistical tests were 2-sided. Baseline quantitative and qualitative variables were analyzed with Wilcoxon-Mann-Whitney U test and Fisher exact test, respectively. Genetic associations were tested by logistic regression with an additive model and adjusted for appropriate clinical covariates unless indicated otherwise, e.g. where logistic regression was not possible. Covariates for logistic regression were derived from each study cohort and represented relevant baseline (clinical and demographic) differences between cases and controls. Covariates for the European Canadian patients (stage 1) analysis included age at start of treatment, cumulative anthracycline dose, tumor type (acute lymphoblastic leukemia, Ewing’s sarcoma and rhabdomyosarcoma) and cardiac radiation therapy. Cumulative anthracycline dose was included as a covariate for the European Dutch patient (stage 2) analysis. Covariates for the combined European patient (stages 1 and 2) and overall combined (stages 1–3) analyses were the same as for the Stage 1 analysis. All effect sizes (odds ratios) were calculated for the minor allele for each SNP.

We performed genetic fine mapping analysis by genotype imputation of additional SNPs not present on the GWAS genotyping platform using BEAGLE 4 with LD and haplotype information from CEU 1000 Genomes population as the reference population. We imputed 1,005,286 additional variants on the chromosome 12 region containing RARG. The BEAGLE allelic R2 was calculated for all imputed SNPs. We then examined evidence of additional genetic associations with ACT in this region based on SNPs with imputed BEAGLE allelic R2 ≥ 0.5, using logistic regression adjusted for stage 1 covariates. LD analyses (r2 and D′) were conducted using the 1000 Genomes CEU reference population.

To control for Type I error, we implemented multiple testing corrections for all genetic association analyses. A threshold of P<1×10−5, indicative of putative genetic associations[47] was implemented for the stage 1 analysis. For the stage 2 replication analysis in European Dutch patients, a Bonferroni-adjusted threshold of P < 0.006 (replication testing of SNPs from 9 LD blocks) was implemented. Further stage 3 replication in independent non-European populations was tested using a threshold of P < 0.05 (1 replicated SNP). Potentially confounding clinical risk factors were identified at P < 0.05.

To prevent spurious genetic associations, all patients within each study population shared the same genetic ancestry, which was self-reported and ascertained by PCA[58]. Then for each study cohort, we computed the genomic inflation factor (λGC) to identify the presence of intra-ethnic fine-scale population structures or admixtures or inflation of the test statistics due to population stratification[59]: Stage 1 cohort (based on 657,694 SNPs) – λGC = 1.021; stage 2 cohort (based on 4516 SNPs) – λGC = 1.014; and stage 3 cohort (based upon 7798 mutual SNPs available to all patients) – λGC = 0.941. We performed meta-analyses of all European patients (stages 1–2) and of all study populations (stages 1–3) using SVS/HelixTree 8.1.1 and calculated the heterogeneity by Cochran’s Q statistics to assess the diversity across the different study populations using Comprehensive Meta Analysis software. The Manhattan plot of −log10 P values and the quantile-quantile distribution were generated using SVS/HelixTree 8.1.1. The regional association plot for the associated genomic region was created using LocusZoom. Linkage disequilibrium plots were created using Haploview and color coded as follows: white (D′ < 1, LOD < 2); blue (D′ = 1, LOD < 2); pink shading (D′ < 1, LOD ≥ 2); bright red (D′ = 1, LOD ≥ 2).

Statistical analyses of the functional data provided in Figure 2 showed that all data (Fig. 2a–c) were normally distributed when assessed by the D’Agostino & Pearson omnibus normality test. Variation coefficients were estimated as: Fig. 2a – RARG WT (23.77%), RARGS427L (31.28%); Fig. 2b – untransfected (9.30%), RARG WT (9.37%), RARG WT + ATRA (7.33%); Fig. 2c – RARG WT (16.48%), RARGS427L (22.69%). F tests showed the variances between groups significantly differed in Fig. 2b and 2c.