Possible modification of BRSK1 on the risk of alkylating chemotherapy-related reduced ovarian function.

STUDY QUESTION: Do genetic variations in the DNA damage response pathway modify the adverse effect of alkylating agents on ovarian function in female childhood cancer survivors (CCS)? SUMMARY ANSWER: Female CCS carrying a common BR serine/threonine kinase 1 (BRSK1) gene variant appear to be at 2.5-fold increased odds of reduced ovarian function after treatment with high doses of alkylating chemotherapy. WHAT IS KNOWN ALREADY: Female CCS show large inter-individual variability in the impact of DNA-damaging alkylating chemotherapy, given as treatment of childhood cancer, on adult ovarian function. Genetic variants in DNA repair genes affecting ovarian function might explain this variability. STUDY DESIGN, SIZE, DURATION: CCS for the discovery cohort were identified from the Dutch Childhood Oncology Group (DCOG) LATER VEVO-study, a multi-centre retrospective cohort study evaluating fertility, ovarian reserve and risk of premature menopause among adult female 5-year survivors of childhood cancer. Female 5-year CCS, diagnosed with cancer and treated with chemotherapy before the age of 25 years, and aged 18 years or older at time of study were enrolled in the current study. Results from the discovery Dutch DCOG-LATER VEVO cohort (n = 285) were validated in the pan-European PanCareLIFE (n = 465) and the USA-based St. Jude Lifetime Cohort (n = 391). PARTICIPANTS/MATERIALS, SETTING,
METHODS: To evaluate ovarian function, anti-Müllerian hormone (AMH) levels were assessed in both the discovery cohort and the replication cohorts. Using additive genetic models in linear and logistic regression, five genetic variants involved in DNA damage response were analysed in relation to cyclophosphamide equivalent dose (CED) score and their impact on ovarian function. Results were then examined using fixed-effect meta-analysis. MAIN RESULTS AND THE ROLE OF CHANCE: Meta-analysis across the three independent cohorts showed a significant interaction effect (P = 3.0 × 10-4) between rs11668344 of BRSK1 (allele frequency = 0.34) among CCS treated with high-dose alkylating agents (CED score ≥8000 mg/m2), resulting in a 2.5-fold increased odds of a reduced ovarian function (lowest AMH tertile) for CCS carrying one G allele compared to CCS without this allele (odds ratio genotype AA: 2.01 vs AG: 5.00). LIMITATIONS, REASONS FOR CAUTION: While low AMH levels can also identify poor responders in assisted reproductive technology, it needs to be emphasized that AMH remains a surrogate marker of ovarian function. WIDER IMPLICATIONS OF THE
FINDINGS: Further research, validating our findings and identifying additional risk-contributing genetic variants, may enable individualized counselling regarding treatment-related risks and necessity of fertility preservation procedures in girls with cancer. STUDY FUNDING/COMPETING INTEREST(S): This work was supported by the PanCareLIFE project that has received funding from the European Union's Seventh Framework Programme for research, technological development and demonstration under grant agreement no 602030. In addition, the DCOG-LATER VEVO study was funded by the Dutch Cancer Society (Grant no. VU 2006-3622) and by the Children Cancer Free Foundation (Project no. 20) and the St Jude Lifetime cohort study by NCI U01 CA195547. The authors declare no competing interests. TRIAL REGISTRATION NUMBER: N/A.

Introduction

Advances in childhood cancer treatment have increased cancer survival rates leading to a growing population of childhood cancer survivors (CCS) (Trama ). Abdominal-pelvic radiotherapy and alkylating agents may compromise ovarian function (Green ; Overbeek ; van der Kooi ) and reduce survivors’ reproductive window. This may manifest as sub- or infertility (Chow ; Anderson ) and a higher risk of premature menopause (Levine ), which in turn may impair quality of life (Langeveld ; van den Berg ; Duffy and Allen, 2009; Carter ; Zebrack ; van der Kooi ). Substantial inter-individual variability in the impact of treatment on ovarian function in similarly treated CCS suggests a role for genetic factors in modifying the association between treatment and the risk of ovarian impairment.

Large-scale genome wide association studies (GWAS) in the general population have identified single-nucleotide polymorphisms (SNPs) associated with age at natural menopause or premature ovarian insufficiency (POI) (Perry ; Stolk et al., 2009; He ; Perry et al., 2013; Day , 2017). These SNPs include variants associated with the DNA damage response (Perry ). Alkylating agents, common chemotherapeutic agents used in childhood cancer treatment, induce apoptosis of cancer cells by damaging DNA and inhibiting cellular metabolisms, DNA replication and transcription (Guainazzi and Schärer, 2010; Kondo ; Fu ). We hypothesized that girls and young women with less efficient DNA damage response systems are more vulnerable to the adverse effects of alkylating agents leading to ovarian dysfunction later in life compared to women with a fully efficient DNA damage repair system.

Serum levels of anti-Müllerian hormone (AMH), produced by the granulosa cells of small growing follicles in the ovaries, are related to age at onset of menopause in healthy women (van Disseldorp ) and can detect ovarian dysfunction prior to both detectible changes in FSH/LH or oestrogen and clinical manifestations of menopause (van Beek ; Nelson ; Anderson ; Dewailly ). In addition, AMH has been demonstrated as a useful and early surrogate marker of reduced ovarian function in cancer survivors (van Beek ; Lie ; Charpentier ; Lunsford ; van den Berg ; van der Kooi ).

Identifying genetic risk factors for treatment-related reduced ovarian function may have clinical implications for risk assessment and medical decision-making regarding fertility preservation in newly diagnosed girls with cancer (van den Heuvel-Eibrink ). The aim of the current study was, therefore, to evaluate whether SNPs in the DNA damage response pathway modify the adverse effect of alkylating agents on ovarian function in CCS.

Materials and methods

Study participants—discovery cohort

CCS for the discovery cohort were identified from the Dutch Childhood Oncology Group (DCOG) LATER VEVO-study, a multi-centre retrospective cohort study evaluating fertility, ovarian reserve and risk of premature menopause among adult female 5-year survivors of childhood cancer (Overbeek ). Data on prior cancer diagnoses and treatments were collected from medical files and information on use of hormones (contraceptives or hormonal replacement therapy) and menopausal status at time of study was obtained from the DCOG LATER VEVO-study questionnaire (Overbeek ). The study was approved by the Medical Ethics Review Committee (IRB protocol number 2006/249, VUmc) and written informed consent was obtained from all participants.

Inclusion and exclusion criteria

Female 5-year CCS, diagnosed with cancer and treated with chemotherapy before the age of 25 years, and aged 18 years or older at time of study were enrolled in the current study. Eligible participants provided a blood sample to quantify AMH levels and extract DNA. Some types of treatment are known to have an invariably extremely detrimental effect on ovarian function. Effects can be so absolute, that this leaves little room for inter-individual variance of the chosen phenotype, as a result of genetic susceptibility. To maximize the potential to detect a role of genetic variation, we excluded survivors who received treatments associated with extensive gonadal toxicity including allogeneic stem cell transplantation, total body irradiation, bilateral ovary-exposing radiotherapy, cranial and/or craniospinal radiotherapy, or bilateral oophorectomy.

Study participants—replication cohorts

PanCareLIFE cohort

PanCareLIFE is a pan-European research project including 28 institutions from 13 countries addressing ototoxicity, fertility and quality of life (Byrne ). This cohort included all adult 5-year female survivors from the PanCareLIFE cohort who were treated for cancer before the age of 25 years and fulfilled all inclusion criteria of this study (van der Kooi ). Demographic, disease- and treatment-related data were collected from medical record files. Approval was obtained from all relevant local review boards and written informed consent from all participants.

St. Jude lifetime cohort

The St. Jude Lifetime Cohort Study (SJLIFE) is a cohort study among 10-year CCS in North America coordinated by the St. Jude Children’s Research Hospital (Memphis, TN, USA) combining treatment data, patient-reported outcomes and clinical assessment (Hudson ). Participants in SJLIFE who fulfilled the inclusion criteria and had blood samples available for AMH and DNA analysis comprised the second replication cohort. Sex hormone use at time of study was documented.

Outcome and outcome definition

The outcome of this study was ovarian function, primarily determined by serum levels of AMH. AMH levels of all three cohorts were determined in the endocrine laboratory of the Free University (VU) Medical Center Amsterdam by an ultra-sensitive Elecsys AMH assay (Roche Diagnostics GmbH, Mannheim, Germany) with an intra-assay coefficient of variation of 0.5–1.8%, a limit of detection (LoD) of 0.01 µg/l, and a limit of quantitation (LoQ) of 0.03 µg/l (Gassner and Jung, 2014).

To account for age-dependency of AMH, participating women in each cohort were divided into four age categories: ≥18–25; ≥25–32; ≥32–40; ≥40 years. These age cut-offs were chosen based on patient numbers, driven by power among the groups, as well as clinical relevance. In each cohort and for each age category, AMH was divided into tertiles with exception of the last age category in which AMH levels varied too little to adequately define tertiles. CCS with an AMH level in the lowest tertile for their age category were defined as having a reduced ovarian function (case), while those with an AMH-value in the highest tertile for their age category were assumed not to have a reduced ovarian function (control). Women over 40 years of age were not considered a ‘case’ based on having an AMH-value in the lowest tertile, but on whether or not they had reported a premature menopause (absence of menses for >12 months before the age of 40) at time of study. No ‘control’ subjects were defined in this age group due to the inability to identify with sufficient certainty those without a reduced ovarian function.

Candidate gene variant selection

SNPs were selected based on a literature search of recently published GWAS that identified loci associated with age at natural menopause (Stolk ; He ; Perry ; van Dorp ). Five GWAS hits in DNA damage response pathways, specifically in the inter-strand cross-link repair pathway, were selected based on the lowest P-value in the largest available GWAS meta-analysis, with the hypothesis that polymorphisms in these regions may increase the gonadotoxic effect of alkylating agents. The selected polymorphisms were in UIMC1 (rs365132), FANCI (rs1054875), RAD51 (rs9796), BRSK1 (rs11668344) and MCM8 (rs16991615). Details concerning the genotype data and quality control protocol are provided in the Supplementary materials and methods file, sections ‘Quality protocol’ and ‘Linkage disequilibrium’.

Alkylating agents

For each survivor, the administered cumulative dose of alkylating agents was quantified using the validated cyclophosphamide equivalent dose (CED) score (Green ). To evaluate the effects of no, low-, medium- and high-dose alkylating agent exposure, the CED score was divided into four categories (0; >0–4000 mg/m2; ≥4000–8000 mg/m2; ≥8000 mg/m2) (Green ). Details on the administered chemotherapeutics, CED score in categories and a fractional polynomial selection procedure for CED score are further discussed in the Supplementary Tables SI, SII, SIII, SIV and SV.

Statistical analyses

Additive genetic associations, with AMH levels based on imputed allelic dosage, were evaluated by logistic and linear regression analyses based on two models: (i) a main effect model; and (ii) an interaction model. Both models evaluated the association between reduced ovarian function and selected SNPs, adjusted for: ancestry and cohort effects using principle components, CED score (four categories using CED of zero as the reference category) (Green ), use of sex hormones (replacement or contraception) at time of study (yes/no), age at time of study (linear regression analysis only) and imputed numbers (0–2) of the alternative allele of the investigated variant (additive effects). The interaction model additionally included an interaction term (SNP*CED category) for genetic variant and CED score categories to evaluate the modifying effect of the variant on the impact of CED score on low AMH levels. Results of linear and logistic regression analyses are presented as regression coefficients (beta) with SE and odds ratios (ORs) with a 95% CI. For linear regression, AMH-levels were log-transformed to adjust for the skewed residuals distribution. Sensitivity analyses performed to assess the robustness of our findings, choices of the model and linkage disequilibrium (Ward and Kellis, 2012) are shown in Supplementary Table SVI. SNPs that showed an association with log-transformed AMH levels or reduced ovarian function in either model, or an interaction effect with CED (P-values < 0.05) were selected for replication of both models. These analyses were conducted using SPSS (Statistical Package for Social Sciences (SPSS) version 24.0.0.1).

Replication and meta-analysis

Findings from the discovery cohort were evaluated in both replication cohorts using identical models, except for sex hormone use at time of study, which was only available in SJLIFE. Data of the discovery and replication cohorts were combined and examined using meta-analytic approaches, in R version 3.5.1, package ‘rmeta’ (R Development Core Team, 2014), the overall P-values for interaction were meta-analysed using Fisher’s method. Pooled estimates based on fixed-effects meta-analysis are presented. In the meta-analysis, P-values <0.01 (0.05/5 gene variants, correcting for multiple testing) were considered statistically significant. Finally, we calculated the cumulative ORs for every genotype per CED category based on the prevalence of a reduced ovarian function for every genotype and every CED category compared to the prevalence of a reduced ovarian function for survivors with a AA genotype treated without alkylating agents, to allow interpretation of the findings.

Results

Discovery cohort

In total, 285 CCS from the DCOG LATER-VEVO cohort participated in the current study (Table I). AMH levels per age category are depicted in Table II. Allele frequencies of the investigated SNPs are depicted in Table III. All SNPs were in Hardy–Weinberg equilibrium (significance level <1*10−7). Results from logistic regression analyses showed a negative association between BRSK1 (rs11668344) and reduced ovarian function (OR 0.56, 95% CI 0.35–0.90; P-value = 0.016) in the main effect-model. In addition, a non-significantly modifying effect of BRSK1 (rs11668344, minor allele frequency 0.34) on the effect of CED ≥8000 mg/m2 on reduced ovarian function (OR 5.02, 95% CI 0.76–33.08; P-value = 0.09) (Table III) was observed in the interaction model. A significant modifying effect of a polymorphism in FANCI (rs1054875) on the effect of CED in the category >0–4000 mg/m2 (OR 9.93, 95% CI 2.35–41.98; P-value = 0.002) was also observed (Table III). Sensitivity analyses of the main analysis did not change the results (Supplementary Tables SVI and SVII). Linear regression analysis showed a significant main effect of the BRSK1 gene variant, but not of the other variants (Supplementary Tables SVIII and SIX). The two SNPs within the BRSK1 and FANCI genes were assessed for replication in the two replication cohorts.

Table I

Characteristics of participating CCS in the discovery and two replication cohorts.

	Discovery DCOG LATER-VEVO (n = 285)	Replication PanCareLIFE (n = 465)	Replication St. Jude Lifetime (n = 391)
Age at time of study (years)
Median (range)	26.1 (18.3–52.4)	25.7 (18.0–45.0)	31.3 (19.1–59.5)
Age at diagnosis (years)
Median (range)	5.8 (0.3–17.8)	10.4 (0.0–25.0)	6.9 (0.0–22.7)
18–25 years	0 (0)	21 (4.5)	16 (4.1)
Time since diagnosis (years)
Median (range)	19.7 (6.7–41.4)	17.0 (5.0–39.1)	23.7 (11.0–46.2)
Diagnosis
Leukaemia	112 (39.3)	109 (23.4)	121 (30.9)
Lymphoma	49 (17.2)	154 (33.1)	70 (17.9)
Renal tumors	37 (13.0)	35 (7.5)	27 (6.9)
CNS tumors	3 (1.1)	12 (2.6)	28 (7.2)
Soft tissue sarcoma	23 (8.1)	31 (6.7)	28 (7.2)
Bone tumors	26 (9.1)	45 (9.7)	34 (8.7)
Neuroblastoma	11 (3.9)	35 (7.4)	36 (9.2)
Other	24 (8.4)	44 (9.6)	47 (12.0)
Radiotherapy
No	251 (88.1)	297 (63.9)	268 (68.5)
Yesa	34 (11.9)	170 (36.1)	123 (31.5)
Thorax	22 (7.7)	88 (18.9)	71 (18.2)
Abdomen (above pelvic crest)	3 (1.1)	12 (2.6)	30 (7.7)
Unilateral ovarianb	0 (0)	9 (1.9)	3 (0.8)
Other	20 (7.0)	61 (13.1)	51 (13.0)
CED score
0	106 (37.2)	161 (34.6)	198 (50.6)
>0–4000 mg/m2	80 (28.1)	103 (22.2)	21 (5.4)
≥4000–8000 mg/m2	52 (18.2)	68 (14.9)	78 (19.9)
≥8000 mg/m2	47 (16.5)	133 (28.6)	94 (24.0)
Hormone use at serum sampling
No	199 (69.9)	232 (49.9)	263 (67.3)
Yes	86 (30.1)	116 (24.9)	128 (32.7)
Oral contraceptive-free day 7	70 (24.6)	3 (0.6)	NA
Anytime during oral contraceptive	NA	94 (20.2)	NA
HRT stop 7	2 (0.7)	20 (4.3)	NA
Anytime, with intrauterine device	14 (4.9)	NA	NA
Unknown	0 (0)	117 (25.2)	0 (0)
Unilateral ovarian oophorectomy
No	284 (99.6)	463 (99.6)	391 (100.0)
Yes	1 (0.4)	2 (0.4)	0 (0)
AMH level
Median (range)	2.5 (<0.01–13.1)	2.1 (<0.01–18.5)	1.8 (<0.01–11.9)
Premature menopause (before age 40) and aged ≥40 years at study,	2 (0.7)	NA	4 (1.0)

Values are represented as the number (%) of women, unless indicated otherwise.

Not mutually exclusive.

Likely in radiotherapy field.

AMH, anti-Müllerian hormone in µg/l; CCS, childhood cancer survivors; CED, cyclophosphamide equivalent dose; CNS, central nervous system; DCOG LATER-VEVO, Dutch Childhood Oncology Group (DCOG) LATER VEVO cohort; HRT, hormonal replacement therapy; NA, not available; PanCareLIFE, PanCareLIFE cohort; St. Jude Lifetime, St. Jude Lifetime Cohort.

Table II

AMH levels in tertiles by age categories.

	VEVO	PanCareLIFE	St. Jude Lifetime
Age 18–25	n = 118	n = 209	n = 72
Lowest AMH tertile	1.08 (0.21–2.14)	0.66 (0.01–1.79)	1.48 (0.15–2.20)
Middle AMH tertile	3.07 (2.16–4.08)	2.51 (1.83–3.39)	2.79 (2.22–3.56)
Highest AMH tertile	5.37 (4.23–13.14)	4.98 (3.41–18.50)	4.91 (3.65–11.90)
Age ≥ 25–32	n = 102	n = 156	n = 143
Lowest AMH tertile	1.32 (0.01–2.14)	0.72 (0.01–1.49)	1.16 (0.01–1.84)
Middle AMH tertile	3.09 (2.15–4.59)	2.33 (1.52–3.26)	2.57 (1.98–3.57)
Highest AMH tertile	6.08 (4.65–12.76)	4.32 (3.27–9.08)	4.87 (3.58–10.48)
Age ≥ 32–40	n = 48	n = 89	n = 107
Lowest AMH tertile	0.36 (0.01–0.80)	0.05 (0.01–0.50)	0.51 (0.01–1.04)
Middle AMH tertile	1.33 (0.91–2.16)	1.19 (0.53–1.90)	1.69 (1.05–2.10)
Highest AMH tertile	3.65 (2.19–9.44)	3.42 (1.93–13.50)	3.27 (2.14–7.70)
Age ≥ 40	n = 17	n = 11	n = 69
No tertiles	0.16 (0.01–1.85)	0.47 (0.01–8.89)	0.09 (0.01–8.73)

Values are represented as the median (minimum–maximum), unless indicated otherwise.

VEVO, DCOG-LATER VEVO cohort.

Table III

Association of single nucleotide polymorphisms with reduced ovarian function and CED-score in DCOG LATER-VEVO discovery cohort.

Gene	Variant	Chrom	Ref.	Alt.	MAF	Model	Variant, interaction term	OR (95% CI)	P-value
BRSK1	rs11668344	19	A	G	0.34	1	rs11668344	0.56 (0.35–0.90)	0.016
							CED: 0	1 (ref)	0.001
							‒ >0–4000	1.43 (0.65–3.11)	0.374
							‒ ≥4000–8000	4.74 (1.92–11.71)	0.001
							‒ ≥8000	5.04 (1.66–15.30)	0.004
							Hormones	2.02 (1.00–4.07)	0.049
						2	rs11668344	0.57 (0.25–1.31)	0.186
							CED: 0	1 (ref)	0.133
							‒ >0–4000	1.94 (0.62–6.07)	0.253
							‒ ≥4000–8000	5.46 (1.32–22.66)	0.019
							‒ ≥8000	1.91 (0.44–8.29)	0.386
							SNP*CED: 0	1 (ref)	0.218
							‒ >0–4000	0.66 (0.21–2.13)	0.489
							‒ ≥4000–8000	0.85 (0.23–3.18)	0.807
							‒ ≥8000	5.02 (0.76–33.08)	0.094
							Hormones	2.01 (0.98–4.14)	0.058
FANCI	rs1054875	15	A	T	0.36	1	rs1054875	1.01 (0.61–1.67)	0.975
							CED: 0	1 (ref)	0.001
							‒ >0–4000	1.37 (0.63–2.95)	0.425
							‒ ≥4000–8000	4.17 (1.73–10.05)	0.001
							‒ ≥8000	4.98 (1.66–14.91)	0.004
							Hormones	1.79 (0.91–3.54)	0.094
						2	rs1054875	0.31 (0.11–0.90)	0.032
							CED: 0	1 (ref)	0.009
							‒ >0–4000	0.32 (0.10–1.06)	0.063
							‒ ≥4000–8000	2.19 (0.60–7.95)	0.235
							‒ ≥8000	3.71 (0.84–16.38)	0.084
							SNP*CED: 0	1 (ref)	0.016
							‒ >0–4000	9.93 (2.35–41.98)	0.002
							‒ ≥4000–8000	3.49 (0.78–15.57)	0.102
							‒ ≥8000	2.00 (0.38–10.44)	0.413
							Hormones	1.83 (0.90–3.73)	0.095
MCM8	rs16991615	20	G	A	0.08	1	rs16991615	0.90 (0.38–2.15)	0.817
							CED: 0	1 (ref)	0.001
							‒ >0–4000	1.37 (0.64–2.94)	0.420
							‒ ≥4000–8000	4.16 (1.74–9.97)	0.001
							‒ ≥8000	4.96 (1.65–14.87)	0.004
							Hormones	1.80 (0.91–3.56)	0.089
						2	rs16991615	0.85 (0.21–3.39)	0.820
							CED: 0	1 (ref)	0.005
							‒ >0–4000	1.36 (0.59–3.14)	0.473
							‒ ≥4000–8000	4.48 (1.73–11.58)	0.002
							‒ ≥8000	3.82 (1.22–11.95)	0.021
							SNP*CED: 0	1 (ref)	0.973
							‒ >0–4000	1.07 (0.14–8.06)	0.950
							‒ ≥4000–8000	0.61 (0.05–6.74)	0.683
							‒ ≥8000	NA	NA
							Hormones	1.89 (0.95–3.75)	0.069
UIMC1	rs365132	5	G	T	0.5	1	rs365132	1.09 (0.70–1.69)	0.720
							CED: 0	1 (ref)	0.001
							‒ >0–4000	1.35 (0.63–2.91)	0.443
							‒ ≥4000–8000	4.18 (1.75–10.00)	0.001
							‒ ≥8000	5.03 (1.68–15.11)	0.004
							Hormones	1.80 (0.91–3.54)	0.090
						2	rs365132	0.79 (0.39–1.61)	0.518
							CED: 0	1 (ref)	0.017
							‒ >0–4000	0.44 (0.11–1.82)	0.257
							‒ ≥4000–8000	4.05 (1.01–16.19)	0.048
							‒ ≥8000	4.83 (0.78–29.90)	0.091
							SNP*CED: 0	1 (ref)	0.265
							‒ >0–4000	2.89 (0.93–8.98)	0.067
							‒ ≥4000–8000	1.04 (0.32–3.39)	0.948
							‒ ≥8000	1.01 (0.17–5.98)	0.988
							Hormones	1.78 (0.89–3.57)	0.104
RAD51	rs9796	15	A	T	0.42	1	rs9796	0.94 (0.62–1.44)	0.787
							CED: 0	1 (ref)	0.001
							‒ >0–4000	1.37 (0.64–2.94)	0.419
							‒ ≥4000–8000	4.17 (1.74–9.99)	0.001
							‒ ≥8000	4.98 (1.66–14.92)	0.004
							Hormones	1.79 (0.91–3.53)	0.092
						2	rs9796	0.92 (0.43–1.97)	0.838
							CED: 0	1 (ref)	0.167
							‒ >0–4000	1.66 (0.52–5.33)	0.397
							‒ ≥4000–8000	4.33 (1.18–15.91)	0.027
							‒ ≥8000	2.34 (0.48–11.42)	0.291
							SNP*CED: 0	1 (ref)	0.546
							‒ >0–4000	0.81 (0.28–2.33)	0.692
							‒ ≥4000–8000	0.94 (0.29–3.16)	0.938
							‒ ≥8000	2.82 (0.52–15.37)	0.230
							Hormones	1.70 (0.85–3.39)	0.135

Alt, alternative allele; Chrom., chromosome; MAF, minor allele frequency; NA, not available; OR, odds ratio; Ref, reference allele; SNP, single-nucleotide polymorphism.

Position based on position build 37 on https://www.ncbi.nlm.nih.gov/snp/. Alt is reported as 0/1/2 (recalculated for presentation only, based on allelic dosage) for CCS with and without reduced ovarian function (see Methods section for details). Model 1: adjusted for principal components, use of hormone use and CED-categories. Model 2: additional to Model 1 interaction term of variant*CED category.

The PanCareLIFE and SJLIFE replication cohorts included 465 and 391 female CCS, respectively (Table I). Consistency of AMH across the three cohorts is depicted in Table II. Table IV shows the combined analysis of both replication cohorts and the final meta-analysis including all three cohorts. Separate findings of the replication cohorts can be found in Supplementary Tables SX and SXI. Full details of the meta-analysis and its heterogeneity are described in Supplementary Tables SXII and SXIII, The overall P-value for interaction between rs11668344 (BRSK1) and CED was 0.018. All three single-cohort analyses suggest a consistent modifying effect for the G allele of rs11668344 (BRSK1) on the effect of CED ≥8000 mg/m2 on reduced ovarian function, although the relatively small-sized discovery cohort did not reach significance for this association. The fixed-effects meta-analysis showed an interaction effect of carrying the G allele of rs11668344 in BRSK1 and an exposure to alkylating agents equivalent to a CED score ≥8000 mg/m2 of 3.81 (95% CI 1.85–7.86, P = 3.0 × 10−4), indicating that the odds of reduced ovarian function increased with an increasing number of G alleles and CED score ≥ 8000 mg/m2.Table V shows the ORs for any genotype per CED category compared to female CCS with the AA genotype and treated without alkylating agents. Female CCS who received alkylating agents equivalent to a CED score ≥8000 mg/m2 had a 2.5-fold higher odds of having an AMH serum level in the lowest tertile with one instead of none G allele of rs11668344 in BRSK1 (genotype AG 5.00 (95% CI 3.27–7.63): AA 2.01 (95% CI 1.31—3.08)) and a 3-fold increased odds with the genotype GG (OR 6.53 95% CI 2.36–18.05).

Table IV

Association of single-nucleotide polymorphisms with reduced ovarian function and chemotherapy in the meta-analyses.

					Replication (PCL+SJLIFE) meta-analysis			Discovery + Replication (VEVO + PCL + SJLIFE) meta-analysis
Gene	Variant	Ref>Alt	Model	variant, interaction	OR (95% CI)	Direction	P-value	OR (95% CI)	Direction	P-value
BRSK1	rs11668344	A>G	2	rs11668344	0.82 (0.54–1.24)	−+	0.349	0.76 (0.53–1.11)	−−+	0.152
				CED: 0	1 (ref)		5.5 × 10−4	1 (ref)		5.6 × 10−4
				‒ >0–4000	0.58 (0.21–1.58)	−−	0.284	0.98 (0.46–2.09)	+−−	0.964
				‒ ≥4000–8000	3.42 (1.52–7.67)	++	2.8 × 10−4	3.83 (1.90–7.74)	+++	1.8 × 10−4
				‒ ≥8000	1.77 (0.18–17.60)	+−	0.627	1.82 (0.40–8.34)	++−	0.442
				SNP*CED: 0	1 (ref)		0.016	1 (ref)		0.018
				‒ >0–4000	3.27 (1.11–9.66)	+−	0.032	1.37 (0.29–6.51)	−+−	0.690
				‒ ≥4000–8000	1.04 (0.44–2.48)	+−	0.922	0.98 (0.48–2.02)	−+−	0.960
				‒ ≥8000	3.63 (1.66–7.95)	++	1.3 × 10−3	3.81 (1.85–7.86)	+++	3.0 × 10−4
FANCI	rs1054875	A>T	2	rs1054875	1.01 (0.65–1.56)	+−	0.977	0.85 (0.57–1.28)	−+−	0.432
				CED: 0	1 (ref)		0.002	1 (ref)		2.0 × 10−4
				‒ >0–4000	0.88 (0.28–2.80)	+−	0.828	0.54 (0.23–1.24)	−+−	0.148
				‒ ≥4000–8000	5.29 (2.08–13.50)	++	4.7 × 10−4	3.91 (1.83–8.33)	+++	4.1 × 10−4
				‒ ≥8000	3.69 (0.37–36.8)	++	0.266	3.70 (0.83–16.6)	+++	0.088
				SNP*CED: 0	1 (ref)		0.869	1 (ref)		0.146
				‒ >0–4000	1.35 (0.46–3.96)	++	0.583	2.76 (1.17–6.53)	+++	0.021
				‒ ≥4000–8000	0.64 (0.29–1.40)	−−	0.264	0.92 (0.46–1.86)	+−−	0.823
				‒ ≥8000	1.03 (0.53–2.03)	++	0.925	1.14 (0.61–2.12)	+++	0.691

PCL, PanCareLIFE cohort; SJLIFE, St. Jude Lifetime Cohort.

Model 2: adjusted for principal components, hormone use (only for VEVO, SJLIFE) and CED−categories and the interaction term of variant*CED category. + = positive association of the SNP with reduced ovarian function in PCL and SJLIFE respectively. − = negative association of the SNP with reduced ovarian function in VEVO, PCL and SJLIFE, respectively.

Table V

OR per genotype of rs11668344 (BRSK1) and CED score on reduced ovarian function, based on prevalence in three cohorts.

	genotype AA		genotype AG		genotype GG
CED in mg/m2	n (%)	OR (95% CI)	n (%)	OR (95% CI)	n (%)	OR (95% CI)
0	51 (40.8)	1 (ref)	36 (40.0)	0.97 (0.63–1.48)	14 (31.8)	0.68 (0.35–1.30)
>0–4000	19 (37.3)	0.86 (0.48–1.53)	19 (38.8)	0.92 (0.51–1.64)	5 (29.4)	0.60 (0.20–1.82)
≥4000–8000	36 (69.2)	3.26 (1.95–5.46)	36 (66.7)	3.48 (2.07–5.87)	7 (43.8)	1.13 (0.41–3.14)
≥8000	43 (58.1)	2.01 (1.31–3.08)	62 (77.5)	5.00 (3.27–7.63)	18 (81.8)	6.53 (2.36–18.05)

n (%) represents the number of cases with reduced ovarian function (% of total) within each genotype group. OR (95% CI) calculated based on the prevalence of a reduced ovarian function for every genotype and every CED category compared to the prevalence of a reduced ovarian function for survivors with a AA genotype treated without alkylating agents.

Linear regression analysis of BRSK1 showed inconsistent associations with AMH in the two replication cohorts, and no significant association was reached in the meta-analysis (Supplementary Table SXIII: beta −0.09, 95% −0.25–0.08). The modifying effect of >0–4000 CED in FANCI (rs1054875) was non-significant in both replication cohorts, and did not reach significance in the meta-analysis (OR 2.76, 95% CI 1.17–6.53, P = 0.02) after correction for multiple testing.

Discussion

This is the first study to assess the influence of genetic factors on alkylating chemotherapy-induced reduced ovarian function, using AMH as a biomarker, and incorporating two independent and identically phenotyped replication cohorts and a meta-analysis. We report a strong modifying effect of a common SNP (minor allele frequency 0.34) in the BRSK1 gene on the toxicity of high dose alkylating agents, resulting in a 2.5-fold increased odds of a reduced ovarian function for CCS carrying one G allele compared to CCS without this allele and a 3-fold increased odds for CCS carrying two G alleles.

One previous single-centre study evaluated the association between ovarian function in CCS with SNPs associated with age at menopause in the general population reporting that the T allele of rs1172822 of the BRSK1 gene was inversely associated with serum AMH levels (van Dorp ). However, this study did not assess interaction between treatment and AMH levels or include validation using replication cohorts. Recently, a SJLIFE GWAS study identified a haplotype associated with an increased risk of premature menopause, especially in the subgroup of CCS who had received pelvic radiotherapy (Brooke ). However, the haplotype is beyond the scope of this study as our population excluded survivors treated with bilateral ovarian radiotherapy due to low inter-individual variation of POI and the haplotype is not associated with DNA damage response genes.

The meta-analysis suggests a strong modifying effect of a G allele of a genetic variant in BRSK1 (rs11668344 A>G) on alkylating agent-related reduced ovarian function. The meta-analysis on reduced ovarian function for the main effect of BRSK1, which is associated with an earlier age at menopause in the general population (Stolk ; He ; Perry ), did not find a significant association as the previous single-centre study reported (van Dorp ). Representing continuous variables such as CED-score in categories may lead to increased type I error for the detection of interaction effects (Royston and Altman, 1994). Supplementary analyses using fractional polynomials (Supplementary Tables SIII, SIV and SV) show that using the available data, estimating more flexible models to potentially avoid these spurious findings, offers inconclusive results due to lack of power, while not contradicting the results found using the pre-defined categories.

Rs11668344 is an intronic variant in THEM150B and an expression quantitative trait locus that alters BRSK1 RNA gene expression in whole blood (P-value = 2.4 × 10−19) (Westra ) and has regulatory histone marks, suggesting a regulatory function. Several mechanisms for the modifying effect of BRSK1 on reduced ovarian function in CCS can be considered. Alkylating agents are known to induce apoptosis of cancer cells by damaging DNA and inhibiting cellular metabolism, DNA replication and DNA transcription (Guainazzi and Schärer, 2010; Kondo ; Fu ). We hypothesize that due to a less efficient DNA damage response system, cancer patients carrying the G allele of rs11668344 in BRSK1 are at an increased risk of the DNA-damaging impact of alkylating agents in healthy tissues most relevant to our outcome studied here, the ovary (Fig. 1). It is plausible that the efficiency of the DNA damage response system becomes crucial upon treatment with alkylating agents amounting to high CED scores.

Figure 1.

Simplified representation of the hypothesized biological plausibility of the effect of DNA damage can be the result of environmental exposure, DNA replication errors but also of chemical exposure. Alkylating agents are known to induce apoptosis of cancer cells by damaging DNA and inhibiting cellular metabolism and DNA replication and transcription (Guainazzi and Schärer, 2010; Kondo ; Fu ). DNA damage response genes (BRSK1 is known to act as a DNA damage checkpoint) have previously been associated with age at natural menopause. Due to a less efficient DNA damage response system, childhood cancer patients carrying the G allele of rs11668344 (BRSK1) may be at an increased risk of the DNA-damaging impact of alkylating agents.

Future research will need to evaluate the relevant expression, which we would expect in granulosa cells or the primordial follicle pool—as opposed to the recruited and selected oocytes that have successfully progressed towards maturation (see also Supplementary file ‘Biological mechanism’).

The identification of this genetic risk factor for alkylating agents-related low AMH levels, if confirmed for other measures of reduced ovarian function, may improve future risk prediction models including more adequate identification of groups with higher or lower risk of chemotherapy-induced ovarian impairment. Upfront fertility preservation programs, including ovarian tissue cryopreservation, would benefit from optimized prediction models as they can be directed to paediatric cancer patients at highest risk for gonadotoxicity for whom the balance of benefits/drawbacks—including ethical considerations—is most beneficial (Warren Andersen, 2018).

A major strength of this study is the inclusion of three independent cohorts which enabled a meta-analysis. As there were some differences between the discovery and the replication cohorts, we performed multiple sensitivity analyses to assess the choices of the model and cohort, which did not change our results. Another strength of this study is the measurement of AMH levels, as a marker for reduced ovarian function, with the same assay at one laboratory, eliminating between-assay differences. Previous studies demonstrated that alkylating agents are strongly associated with risk of reduced ovarian function as measured by decreased AMH levels in female CCS (Anderson ; Thomas-Teinturier ; van der Kooi ; van den Berg ). By using AMH levels as a marker of ovarian function, this study included a fairly substantial number of cases likely at increased risk of reduced fertility or a shorter reproductive window. However, while low AMH levels can also identify poor responders in assisted reproductive technology (Iliodromiti ; van Tilborg ), it needs to be emphasized that AMH remains a surrogate marker of ovarian function. The implications of low AMH on natural fertility and reproductive lifespan are under continuing debate. While in the general population AMH has proven to be a valuable predictor of menopause, apart from age (van Disseldorp ; Tehrani ; Freeman ; Dolleman ; Depmann ), current prediction models have not been designed to predict the extremes of menopausal age (Depmann ). Validation using data collected long-term and using more definite and direct endpoints such as age at menopause, POI, or fecundity is needed to facilitate translation into clinical practice. In addition, larger cohorts would benefit the power of statistical tests.

In conclusion, this study presents data suggesting that high dose alkylating chemotherapy-induced reduced ovarian function in female CCS is strongly modified by a common DNA variant (rs11668344) of the BRSK1 gene. This is the first time a genetic risk factor has been described to modify the effect of chemotherapy on long-term ovarian function in three independent cohorts. This finding may serve as a starting point for further research working towards individualized counselling regarding treatment-related risks and fertility preservation services in children with cancer as well as young adult survivors.

Data availability

The data underlying this article cannot be shared publicly due to ethical reasons and privacy of individuals that participated in the study. The data will be shared on reasonable request to the corresponding author, and after consultation of data and ethics committees of the three separate cohorts.

Click here for additional data file.

Click here for additional data file.

Click here for additional data file.

Click here for additional data file.

Click here for additional data file.

Click here for additional data file.

Click here for additional data file.

Click here for additional data file.

Click here for additional data file.

Click here for additional data file.

Click here for additional data file.

Click here for additional data file.

Click here for additional data file.

Click here for additional data file.

Click here for additional data file.