Identification of genetic variants associated with skeletal muscle function deficit in childhood acute lymphoblastic leukemia survivors

Introduction

Childhood acute lymphoblastic leukemia (ALL) accounts for nearly 25% of all paediatric malignancies.^1,2 Optimization of multi-agent risk-adapted protocols has improved the 5-year survival rate, reaching by 2017 an estimated 90% in children less than 14 years old.³ This therapeutic success introduces new challenges, as 60% of childhood ALL survivors now experience a broad spectrum of late adverse effects (LAE) due to aggressive chemo- and radiation-therapy administered during critical child development stages.^4,5 Childhood ALL survivors are more likely to develop a metabolic syndrome, cardiotoxicity, neuropsychological deficit and musculoskeletal morbidity — a cluster of interconnected factors that globally contribute to severe impairment of their quality of life.^6–11

Clinical evaluation of childhood ALL survivors has revealed skeletal muscle size, density and function deficits as compared with healthy controls.¹² Evidence suggests acute oxidative stress induced alteration of functional muscle mitochondria, motor neurons and satellite cells could contribute to these functional deficits.¹³ Muscle morbidity could lead to long-term mobility and autonomy deficits, such as impaired walking and gait or a diminished capacity to carry out daily activities.^14–16 As lipid oxidation and glucose consumption is primarily performed in skeletal muscle, it is likely that treatment-induced muscle function deficit ultimately contributes to obesity, and cardiovascular or metabolic syndrome-related morbidity in ALL survivors.¹⁷

While molecular mechanisms involved in skeletal muscle atrophy remain unclear, sarcopenia as a response to glucocorticoid (GC) signalling, irrespective of cancer treatment, has been associated with various genetic markers and molecular pathways in the adult population. Evidence shows that GCs act as up-regulators of atrophy-related genes such as the Notch signalling pathway and constituents of the E3 ubiquitin ligase cascade.¹⁸ Likewise, mTOR (mechanistic target of rapamycin kinase) and AMPK (AMP-activated protein kinase) pathways are involved in GC-induced delays of muscle development in an animal model, while GCs increase adiposity in muscle tissue by affecting IL-4 (Interleukin 4) regulation of fibro/adipogenic progenitors.^19,20

Genetic susceptibility to long-term skeletal muscle morbidity in childhood ALL survivors has yet to be documented. To address this issue, we present correlations between relevant clinical outcomes and genotypic profiling obtained through whole exome sequencing (WES) of ALL survivors.

Identification of genetic markers with predictive potential would allow further risk stratification of patients and personalisation of LAE prevention programs.

Methods

Study participants

The study group is composed of 242 survivors from the PETALE cohort, diagnosed with childhood ALL since 1987 at the Sainte-Justine University Hospital Center (SJUHC) in Montreal.¹⁰ Eligibility criteria include an age at diagnosis of 19 years or less, a 2-year adherence to a Dana-Farber Cancer Institute (DFCI)-ALL treatment protocol, Caucasian descent, and a post-remission delay of at least 5 years. Patients who relapsed, received hematopoietic stem cell transplants, suffered from concomitant syndromes, congenital bone disease, sickle-cell anemia or Legg-Calvé-Perthes disease, or were given osteotoxic drugs outside of their ALL treatment were considered non-eligible for the study. All patients fulfilling eligibility criteria and who consented to participate were included in this study.¹⁰ A written informed consent was obtained in accordance with the principles laid down by the Declaration of Helsinki. The SJUHC Institutional Review Board approved the research proposal.

Muscle force and power phenotypic data

Study participants underwent comprehensive clinical, biological and psychosocial investigations carried at the SJUHC and the Montreal Shriners Hospital for Children (MSHC), in the course of a day, between 2013-02-12 and 2016-05-19.

We used the Leonardo Mechanograph® Ground Reaction Force Plate (Novotec Medical GmbH, Pforzheim, Germany) as the measurement device for peak muscle force and power data acquisition. Data were analysed using the Leonardo Mechanography GRFP Research Edition® software, which measures body mass during quiet stance, and calculates the instantaneous vertical acceleration upon jumping, by recording the force as a function of time. The acceleration is time integrated to obtain the instantaneous vertical velocity. Peak power (Pmax) designates the highest value of force times velocity during the initial acceleration phase, while peak force (Fmax) refers to the maximum vertical force applied to the platform. Each participant performed a single two-legged jump and multiple one-legged hopping.²¹ Weight-adjusted Pmax and Fmax data were collected for each participant during a single two-leg jump and transformed into age- and gender-specific z-scores by comparison with a normal population.²² A z-score of −2 was used as a threshold to discriminate between affected and unaffected individuals.

Sequencing and candidate gene selection

A list of candidate genes was generated through a review of the literature. Based on the functional bone-muscle local unit theory, genes involved in both muscle and bone molecular pathways were targeted. Studied pathways include Wnt-β catenin signalling, osteoclastogenesis, osteoblast and osteocyte maturation, Vitamin D signalling, and pathways involved in glucocorticoid and methotrexate (MTX) action.^23–30 Additional candidate genes were selected via a Gene Ontology Enrichment Analysis and Visualisation Tool (Gorilla).³¹ A total of 1039 candidate genes were chosen.

DNA was extracted from whole blood obtained during phenotyping evaluation of PETALE participants. Whole exomes were captured in solution with Agilent’s SureSelect Human All Exon 50Mb kits and sequenced on either Life Technologies SOLiD System 4.0 (mean coverage = 40X) or Illumina HiSeq 2500 platform (mean coverage = 113.1X) at SJUHC integrated clinical genomic centre in paediatrics. Reads were aligned to the hg19 reference genome using SOLiD LifeScope software for the SOLiD samples and BWA-MEM for the Illumina samples.^32,33 PICARD was used to mark PCR duplicates and collect sequencing quality control metrics.³⁴ Germline variant calling was performed using the Haplotype Caller and quality score recalibration was performed using Variant Recalibrator, both implemented in the Genome Analysis Tool Kit (GATK).³⁵ Variants were selected based on the variant quality score (VQSR = PASS) and minimum depth of coverage (DP ≥ 10). WES data were available for 222 participants. Common and rare variants were defined respectively by their minor allele frequency (MAF) in the reference population as ≥5% and <5%, as based on 1000genome and ESP6500 (European American population) datasets.³⁶ A functional annotation from ANNOVAR was used to apply filters to genotype data so that only variants with missense (at least potentially damaging, Polyphen2 scores ≥ 0.85 and Sift scores ≤ 0.1) or nonsense variations, or variations in splicing sites, were further considered.³⁷

Of those, common variants either exceeding a missing rate of 20%, not in the Hardy-Weinberg Equilibrium (HWE, p<0.001), or having a pairwise linkage disequilibrium (LD, r2>0.8) were excluded, leaving 253 common SNPs located in 188 genes. The rare variants were analysed in an exploratory fashion using a collapsing approach whereby minor alleles of at least 2 potentially functional variants within one gene, each in HWE and with missing rate <20%, were analysed together, resulting in an analysis of 655 genes.

Statistical analysis

We used a whole genome data analysis toolset (PLINK^©) to quantify associations between common SNPs and phenotypic datasets applying an allelic ratio test for dichotomous variables and a linear regression for continuous phenotypes. Quantitative datasets were tested for normality through an assessment of the distribution of residuals. For rare variant associations, the SKAT-O test (Optimal Sequence Kernel Association Test) in R Studio© software, which collapses variants within the same gene, was used.³⁸ Resulting P-values were adjusted for multiple testing using a Benjamini-Hochberg correction method. A maximum false discovery rate (FDR) of 5% was chosen as a cut-off value for significance.

Selective genotyping of top-ranking SNPs (FDR<5%) following WES data analysis was carried out for all participants on the Sequenom platform at the McGill University and Génome Québec Innovation Centre. DNA samples from all 242 participants were available for genotyping. Depending on the studied SNP, genotyping was inconclusive for 3–14 participants. Concordance percentages between sequencing and genotyping results for top-ranking signals were recorded. Resulting genotypes and phenotypic data were analysed by IBM SPSS Statistics for Mac Version 24.0 (IBM Corp. Armonk, NY) with the chi-square or the Fisher exact tests, using appropriate genetic models. Rare SNPs were analysed through a collapsing approach with iterative exclusion of a single SNP per iteration, allowing weighting the SNP contribution to association signals.

Logistic regression models (in IBM SPSS) were used to assess an influence of genotypes on studied outcomes in the presence of covariables that significantly correlated with studied phenotypes and included sex, risk status either high (HR) or standard (SR), and body mass index (BMI). BMI was included as a cofactor as adipose tissue accumulation in obese individuals is thought to interfere with muscle function through aberrant adipokines production.³⁹ A participant BMI was classified as normal (BMI ≤ 24.9 kg/m²), overweight (25 kg/m² ≤ BMI < 30 kg/m²), or obese (BMI ≥ 30 kg/m²), in accordance with the adult classification scale.

Due to the DUOX2 gene role in thyroid hormone synthesis, the polymorphism rs2001616 located in DUOX2, for which significant associations with musculoskeletal deficit were found, was analysed for a correlation with TSH and T4 levels available for 235 and 231 PETALE participants, respectively.

To allow better characterisation of the impact of genotypes in relation to Fmax and Pmax, stratified analyses were also performed according to sex, risk group, and BMI at visit (obese vs non-obese).

Results

Cohort characteristics

Fmax measurements were available for 215 participants (Table 1), of whom 40% had a z-score below −1 while 13% had a z-score ≤ −2. As for Pmax, results were available for 226 participants: 55% and 27% with z-scores ≤ −1 and ≤−2, respectively. Thus, the PETALE cohort is more affected in terms of muscle morbidity than the normal population (P<10⁻⁴).

Table 1 Cohort characteristics and distribution of study outcomes

Association study with WES data

Association analyses were performed between functionally predicted common and rare variants derived from selected candidate genes obtained from WES data. Among common SNPs, rs41270041 located in the ADAMTS4 gene, was positively associated with Fmax z-scores ≤ −2 (P=10⁻⁴), whereas rs2001616 in the DUOX2 gene was significantly associated with Pmax, analysed both as a categorical (z-score ≤ −2) and continuous variable (P=10⁻⁵ and P=4.6×10⁻⁵, respectively. Table S1).

Rare-SNP analyses identified top-ranking ALOX15 (P=5.7×10⁻⁶) gene, associated with continuous Pmax z-scores. The signal originated from three collapsed rare variants in the ALOX15 gene (Table S2).

Confirmation by genotyping and association analyses of obtained genotypes

Participants were screened for top-ranking SNPs, including 2 common SNPs (rs41270041 in ADAMTS4 and rs2001616 in DUOX2) and 3 rare variants (rs113604586, rs34210653 and rs11568142) in ALOX15. Concordance between WES and genotyping data ranged from 92% to 95%. Given that similar results were obtained in association analyses for continuous Pmax z-scores and z-scores ≤ −2, only the latter are presented. The rs2001616 in the DUOX2 gene was associated with Pmax according to a dominant model of inheritance (P=2×10⁻⁴, Table 2). Carriers of a minor G allele were more frequent among participants with a z-score ≤ −2 as compared with those carrying the AA genotype (Odds ratio (OR)=3.4, 95% confidence interval (CI)=1.7–6.6). Likewise, the minor C allele of the rs41270041 in ADAMTS4 was associated with Fmax (z-score ≤ −2) according to a dominant model (OR=2.8, 95% CI, 1.2–6.8, P=0.02, Table 2). We next explored whether this association is influenced by clinical and treatment characteristics. Sex, prognostic risk groups, cranial radiation therapy (CRT), BMI, as well as cumulative GC and MTX dose were all significantly associated with Pmax Z score ≤ −2 (P≤0.02, Table S3). Pmax Z score reduction was seen more frequently in female survivors, in obese subjects and participants who were treated according to the high risk protocol. A significant association was noted as well with cumulative MTX and GC doses and with CRT. This reflected more frequent CRT, higher GC doses, but lower MTX doses administered to patients assigned to the high risk group (P<0.00001, Table S4). No significant association with demographics or treatment characteristics was found for Fmax. A logistic regression model was thus derived only for DUOX2 rs2001616, which remained significant in the multivariate model (P=0.003, Table 3) that included, besides genotype, sex, BMI at visit, and risk status (as a surrogate for CRT and higher GC doses). Subsequent stratification analyses were performed and showed that the association of both the DUOX2 rs2001616 and ADAMTS4 rs41270041 polymorphisms with phenotypes was modulated by demographics and treatment characteristics. The association signal for DUOX2 rs2001616 was limited to females (P=0.001), non-obese participants (P=0.0001) and patients assigned to the HR group (P=0.004) (Figure 1), whereas that of the rs41270041 in ADAMTS4 was more apparent among males (P=0.02) and participants that received the HR treatment (P=0.0001) (Figure 2). Given the strong correlation between HR and exposure to CRT (Table S4), as well as with higher GC cumulative doses, the similar modulation of the effect is seen when CRT or GC doses are used as stratification variables (data not shown).

Table 2 Top-ranking common SNPs, genotyping data, and associated phenotypes

Table 3 Regression model

Figure 1 Stratification result of association for DUOX2 rs2001616. Frequencies of AA and AG+GG genotypes (gray bars) and unaffected (black bars) individuals, as defined by Pmax Z≤−2 SD and Z>−2 SD, respectively. P-value with OR (odds radio) and 95% CI (confidence interval) is presented for each comparison. Number of patients in each genotype group is indicated in brackets with percentages. BMI; non obese: <30 kg/m2, obese: ≥30 kg/m2.

Figure 2 Stratification result of association for ADAMTS rs4127041. Frequencies of GG and GC+CC genotypes for affected (grey bars) and unaffected (black bars) individuals, as defined by Fmax Z≤−2 SD and Z>−2SD, respectively. P-value with OR (odds radio) and 95% CI (confidence interval) is presented for each comparison. Number of patients in each genotype group is indicated on each bar.

An analysis of the four different combinations of the rs34210653, rs113604586 and rs11568142 rare variants in the ALOX15 gene revealed that the association signal resulted from the rs113604586 and rs11568142 combination, with a lowest P-value of 0.001 (Table 4). The low frequency of rare variant carriers precluded a stratified analysis.

Table 4 Top-ranking rare SNPs, genotyping data, and associated phenotypes

Discussion

Studies exploring the role of a genetic component in LAE development in childhood ALL survivors are limited, and this study is the first to focus on skeletal muscle function deficit through a genomic approach.^40,41

Our analysis revealed an association of rs2001616 in the DUOX2 gene with skeletal muscle Pmax deficit. Survivors carrying a minor allele had a higher probability of Pmax reduction than other genotype groups. DUOX2 (Dual oxidase 2) is involved in muscle and bone metabolisms through its role in thyroid hormone synthesis. The encoded protein acts as a functional constituent of a peroxide generating system located at the apical membrane of thyroid follicular cells.⁴² Several DUOX2 polymorphisms have been associated with congenital hypothyroïdia, and inactivating mutations of this gene have recently been shown to play a causal role in very early onset inflammatory bowel diseases.^43,44

Thyrotoxicosis is known to be a risk factor for hypercalcemia, osteoporosis and increased bone turnover.⁴⁵ A study investigating correlations between thyroid hormone levels and skeletal muscle regeneration found that high serum T3 and T4 levels could interfere with the inflammatory process following muscle injuries through up-regulation of macrophage accumulation and cell infiltration.⁴⁶ The detected effect of rs2001616 polymorphisms could therefore be due to up-regulation of thyroid hormone production in carriers of the variant allele.

We thus evaluated correlations between the DUOX2 rs2001616 variant and thyroid function in study participants. TSH and free thyroxine (T4) levels were available for 232 and 229 participants with genotype data. A significant association (P=0.02) was found between free T4 values and the rs2001616 polymorphism, carriers of a variant allele having a higher median T4 level (12.8 vs 12.2 μmol/L). Associations with TSH were statistically insignificant. The clinical significance of this result is limited, however, as none of the study participants was diagnosed with dysthyroïdia as based on TSH and T4 levels following their ALL treatment. Furthermore, no significant association was found between thyroid hormone levels and skeletal muscle function phenotypes.

An association with altered skeletal muscle function (Fmax) was also demonstrated with rs41270041, located in the ADAMTS4 gene. ADAMTS4 (ADAM metallopeptidase with thrombospondin type 1 motif 4) is a protease-coding gene responsible for the degradation of aggrecan, a cartilage-specific proteoglycan core protein.⁴⁷ This gene has been shown to contribute to osteoarthritic disorders and tendon injuries.^48,49 An increased susceptibility to tendinopathies or painful joints could lead to limb immobilization and thus alter muscle mass through lower muscle strain. The ADAMTS4 protein has also been shown to be a risk factor for atherosclerosis plaque vulnerability in mice models.⁵⁰

ALOX15 codes for arachidonate 15-lipoxygenase, a pro-inflammatory protein involved in the arachidonic acid metabolism pathway. It has been reported to play a role in inflammatory response and immunity, and in ischemic heart disease.^51,52 Studies carried on mice models suggest that ALOX15 could also affect bone mass through its negative regulation of maximum bone mineral density.⁵³ Associations have been detected between ALOX15 variants and bone micro-architecture in young women, while other variants in this gene seem to predispose to fractures in older women.⁵⁴ We detected a significant association for two rare SNPs included in the ALOX15 locus, justifying further study of this gene in related contexts.

A likely explanation why top-ranking SNPs did not all correlate with both Fmax and Pmax phenotypes lies in the distinct skeletal muscle functions reflected by these two measurements. Muscle force represents the ability to accelerate a given mass (Force = Mass × Acceleration) whereas power is the ability to produce force at a given velocity (Power = Force × Velocity). Previous studies have demonstrated that peak power, but not peak force, depends on muscle fiber type.⁵⁵ On the other hand, evidence suggests a strong correlation between peak muscle force and muscle size as defined by cross-sectional area measurements.^56,57 We can thus hypothesize that the ADAMTS4 gene could modulate molecular pathways involved with skeletal muscle size, while ALOX15 and DUOX2 would affect muscle fiber contractile properties.

Associations of skeletal muscle deficits with common variants were more significant in groups known to suffer more from treatment-related complications (e.g. CRT, HR, and high GC doses groups). In fact, various studies demonstrated that muscle function is substantially more degraded in patients exposed to a highly aggressive treatment regimen, including higher glucocorticoid doses and CRT, which may explain the interaction with risk status noted for associated SNPs of DUOX2 and ADAMTS4 genes.^58,59

A study carried in the St Jude Lifetime Cohort demonstrated severe anterior pituitary deficit related to CRT exposure in childhood cancer survivors.⁶⁰ However, none of the patients in the PETALE cohort were diagnosed with secondary clinical hypopituitarism, which may be due to lower overall CRT doses received. It remains possible that a higher proportion of participants in the CRT-positive group suffered from CRT-induced skeletal muscle deficit due to alternative mechanisms involved with neuroendocrine signalling. Indeed, low dose ionising radiation therapy has been demonstrated to rapidly alter synaptic and mitochondrial signalling pathways in murine models, which could account for poorer outcomes, particularly apparent in patients with a given genetic background.⁶¹

The impact of obesity seems to overcome genetic predisposition, as a higher risk of skeletal muscle function deficit in the case of the DUOX2 SNP was seen only among non-obese patients. Studies investigating muscle pathologies in obese individuals have found that dysregulated adipokines production in hypertrophied adipose tissue decreased the expression of contractile proteins in myotubes, pointing towards an alternative mechanism for poorer skeletal muscle function in this group.³⁹

A limitation of this study results from our representation of skeletal muscle function deficit through muscle power and force phenotypes. Additional measures of muscle quality, architecture, and function and myosteatosis such as muscle cross-sectional area, lean body mass, and metabolic indices could have been used to complement the study phenotypes and better grasp the issue of sarcopenia in childhood leukemia survivors.⁶²

It is not excluded that variables not considered in this study, such as nutritional status, level of physical activity, as well as concomitant endocrine or chronic inflammatory conditions, could have been confounding factors in genetic associations. Some of these factors were not included due to the subjectivity of self-reported phenotypic data (e.g. level of physical activity), and others because none of the PETALE cohort participants presented evidence of such conditions (e.g. hypopituitarism). Likewise, although we controlled the analyses as much as possible for potential major confounders, a number of covariables could have still influenced the analyses, such as wide age range, duration of follow-up, or different complications that occurred during the treatment.^63,64

It is worth noting that the limited sample size might have reduced the accuracy of the stratified analyses, and may have affected the detection of lower risk effects (the lowest OR that can be detected for common variants after correction for multiple testing in WES dataset is 2.6). On the other hand, cohort homogeneity with respect to treatment and ethnic background reduced the impact of confounding variables on the study results.¹⁰ While the candidate gene approach is generally biased towards rational gene selection, the investigation of many such genes allowed the identification of novel variants associated with genetic predictability of skeletal muscle deficit in the targeted population. The analyses of WES data were corrected for multiple testing using a 5% FDR, which in this dataset corresponded well with Bonferroni-adjusted P-values for the number of common variants tested or the number of genes in rare variants testing.

A replication analysis should follow this study to validate findings, whereas future studies conducted at the exome or genome level may detect in an unbiased fashion new genetic predictors underlying skeletal muscle deficit in childhood ALL survivors. Together, they may explain the observed heterogeneity of these late-adverse effects in ALL survivors, and could eventually lead to personalised prevention strategies.