Genetically determined telomere length and multiple myeloma risk and outcome

Telomeres are involved in processes like cellular growth, chromosomal stability, and proper segregation to daughter cells. Telomere length measured in leukocytes (LTL) has been investigated in different cancer types, including multiple myeloma (MM). However, LTL measurement is prone to heterogeneity due to sample handling and study design (retrospective vs. prospective). LTL is genetically determined; genome-wide association studies identified 11 SNPs that, combined in a score, can be used as a genetic instrument to measure LTL and evaluate its association with MM risk. This approach has been already successfully attempted in various cancer types but never in MM. We tested the “teloscore” in 2407 MM patients and 1741 controls from the International Multiple Myeloma rESEarch (IMMeNSE) consortium. We observed an increased risk for longer genetically determined telomere length (gdTL) (OR = 1.69; 95% CI 1.36–2.11; P = 2.97 × 10−6 for highest vs. lowest quintile of the score). Furthermore, in a subset of 1376 MM patients we tested the relationship between the teloscore and MM patients survival, observing a better prognosis for longer gdTL compared with shorter gdTL (HR = 0.93; 95% CI 0.86–0.99; P = 0.049). In conclusion, we report convincing evidence that longer gdTL is a risk marker for MM risk, and that it is potentially involved in increasing MM survival.


Introduction
Multiple myeloma (MM) is a plasma cell malignancy that arises from a single clone of malignant plasma cells (PCs) in the bone marrow. It has a worldwide incidence of 2.1/100,000 new cases every year, ranking as the second most common hematological cancer 1 . MM has usually a late onset with a mean age at diagnosis around 60 years 1 .
There are overwhelming evidences that genetic variability influences the risk of developing MM [2][3][4][5][6][7][8] . In the recent years the importance of genetics has also emerged in MM response to treatment and survival [9][10][11][12][13][14][15] . Alongside polymorphic variants, another emerging marker of susceptibility and prognosis for several diseases is telomere length. Telomeres are specialized structures that cap the end of chromosomes and are involved in several processes like cellular growth, chromosomal stability, and proper segregation to daughter cells 16 . Telomerases are the most important enzymes involved in telomere replication and their dysfunction is considered a cancer hallmark 17 . Telomere length measured in leukocytes (LTL) has been used to investigate the susceptibility to different cancer types including several hematological malignancies [18][19][20][21][22][23][24][25] , among which MM for which only two studies of small size have been attempted, a prospective cohort study and a retrospective case control study 21,26 . In a meta-analysis conducted by Xu and colleagues the role of LTL in relation to survival in several cancer types has been investigated 27 . However, MM was not among the tumors investigated by the authors and the only study conducted on MM survival consisted of five patients where telomere length was measured in bone marrow cells, comparing malignant vs. non-malignant cells 28 .
Telomere length is highly correlated across tissues 29,30 , it is reasonably stable during time, as assessed in a longitudinal study 18 , and therefore it is considered a valid surrogate for the measure of telomere length in specific tissues.
However, LTL measurement is prone to heterogeneity due to sample handling (DNA extraction methods used, storage of blood and/or DNA, elution buffer used) 31 . In addition, the analysis of LTL is sensitive to other confounders, such as the epidemiologic design of the study (retrospective vs. prospective) and, in studies on cancer, whether the subjects have received chemotherapy or not. LTL is genetically determined (gdTL) with an estimated heritability that ranges from 0.44 to 0.70 32,33 . Recently genome-wide association studies (GWAS) have identified 11 SNPs 23,34-36 that, collectively, explain 2.28% of LTL variability and can be used as a genetic instrument in a Mendelian randomization fashion. This approach has been already successfully tested in relation to the risk of developing various cancer types 23,37-45 , but was never done in MM. We have tested the genetic score (henceforth called "teloscore") in the context of the International Multiple Myeloma rESEarch (IMMeNSE) consortium to investigate the relationship between gdTL and MM risk and survival.

Study population
This study was conducted within the IMMeNSE consortium, that has been extensively described elsewhere 46 . We used DNA samples from 2407 MM patients and 1741 controls for whom information on sex, age (age at diagnosis for the cases/age at recruitment for the controls), and country of origin (Denmark, France, Hungary, Israel, Italy, Poland, Portugal, and Spain) was collected. Cases were defined by a confirmed diagnosis of MM according to the International Myeloma Working Group (IMWG) criteria 47 . In addition, clinical data such as disease stage according to Durie-Salmon (DS) staging system 48 and/or International Staging System (ISS) 49 , first-line treatment received, based on bortezomib/immunomodulatory drugs (defined as "new treatments") or any other regimen, such as vincristine/adriamycin/dexamethasone or melphalan/ prednisone (defined as "old treatments") and overall survival (OS) were retrospectively collected from medical records. Information about DS stage was available for 1235 cases, while the stage according to ISS was available for 1022 cases, and for 984 cases information on both systems was available. Controls were selected among the general population as well as among hospitalized subjects or among blood donors in the same geographic areas where the cases were collected. Characteristics of the study population are summarized in Table 1. The
Genotyping DNA samples were extracted from whole blood. Genotyping was conducted in 384-well plates, using TaqMan (ABI, Applied Biosystems, Foster City, CA, USA) technology according to the manufacturer's recommendations. For quality control purposes around 10% of the samples were duplicated in order to check the concordance between the genotypes. Case/control status was unknown to the person performing the genotyping. Fluorescent signals of the genotyping assays were read on PCR plates by a spectrophotometer (Ω, BMG LABTECH, Ortenberg, Germany) and the software KlusterCaller (LGC Group, Teddington, UK) was used to determine genotypes.

Teloscore computation
The "teloscore" was computed as follows: for each SNP the number of alleles associated with longer telomeres (according to the results of the literature) was counted and added up for each study subject, resulting in the unweighted score. Therefore, considering 11 SNPs, the unweighted score could assume any value between 0 (shortest telomeres) and 22 (longest telomeres). We then created a weighted score for each study subject. First, we took from the literature estimates of the per-allele effect on LTL in base pairs for each SNP (Supplementary Table  1). Then, we multiplied for each SNP the number of alleles associated with longer telomeres by the per-allele effect on LTL in base pairs that were then summed up for each study subject. The weighted score thus represents the estimated difference in telomere length, measured in base pairs, attributable to the SNPs under investigation. The concordance between the teloscore built with SNPs and the relative telomere length measured with real time quantitative PCR was described in detail elsewhere 37 . Only a subset of the study subjects had a 100% SNP call rate (1769 cases and 1207 controls). For survival analyses, unweighted and weighted teloscores were built with all the individuals with complete data on stage, therapy received and OS. The number of cases with 100% SNP call rate was 1000 with DS staging system information and 813 with ISS.
In order to be able to compute comparable score values for all study subjects, we also considered average values for each score. Supplementary Table 2 shows examples of how the teloscores were generated.

Statistical analysis
The association between the individual SNPs and MM risk was tested using unconditional logistic regression computing odds ratios (ORs) and 95% confidence intervals (CIs). We used an allelic (log-additive) and a codominant model of inheritance. The threshold for statistical significance was therefore p = 0.05/(11 SNPs x 2 models) = 0.05/22 = 0.0023.
For both teloscores (weighted and unweighted) we calculated quintiles based on the distribution of values in the controls. The association between the teloscores and MM risk was tested with logistic regression computing OR and 95% CIs.
All analyses were adjusted for age, sex, and country of origin.
Survival analysis to assess the effect of the SNPs and the teloscore on OS of MM patients was done using Cox regression, calculating hazard ratio (HR), and 95% CI. OS was defined as the time between MM diagnosis and last follow up for alive patient or death. Survival analyses were carried out with individual SNPs and with both teloscores (weighted and unweighted). We built the teloscores calculating quintiles based on the distribution of values in all the individuals with staging, therapies and OS data. The analyses were adjusted for age, sex, country of origin, MM stage (according to DS or ISS), and first-line therapy.
To better assess the causality of the association between LTL and MM onset and progression we used a Mendelian randomization approach considering the 11 SNPs as genetic instruments of the exposure, and risk of developing MM and MM OS as outcomes of interest. We utilized the inverse variance weighted (IVW) regression because on the premise of absence of directional pleiotropy, it provides consistently robust causal effect estimates 50,51 . The Egger (MR-Egger) estimation method was used to test for possible pleiotropic effects of the genetic instruments 50,51 . Finally to assess evidence of heterogeneity between the polymorphisms, the Cochran Q′ heterogeneity and Bowden I 2 GX statistics were computed to obtain an estimate of instrument strength for the MR-Egger analysis 52 . Additional information has been given elsewhere 50,51 .

Bioinformatic tools
Several bioinformatic tools were used in order to identify functional relevance for the polymorphic variants showing the most promising associations with risk of developing MM. In particular, RegulomeDB (https:// www.regulomedb.org/regulome-search/) 53 was used to identify the regulatory potential of the region nearby each SNP while the GTEx portal web site (https://www. gtexportal.org) was used to identify potential associations between the SNP and expression levels of nearby genes (eQTL) 54 .

Results
Characteristics of the study population are summarized in Table 1 We observed several associations between individual SNPs and MM risk. In particular, two polymorphisms showed a statistically significant association considering the Bonferroni-corrected threshold (p < 0.0023). TERC-rs10936599-T was associated with decreased risk of developing MM (OR homozygous = 0.55; 95% CI 0.39-0.76; P = 2.87 × 10 −4 ), while OBFC1-rs9420907-C was associated with increased MM risk (OR heterozygous = 1.32; 95%CI 1.12-1.55; P = 7.78 × 10 −4 ). All the results of the risk analysis considering the SNPs individually are shown in Table 2.
Association of the "teloscore" with MM risk Supplementary Table 2 illustrates how the teloscores were generated. We tested the association between the unweighted and weighted teloscore with MM risk in individuals that had 100% call rate only (1769 cases and 1207 controls). We observed an increased risk for individuals that had longer gdTL. The association was similar in terms of estimates and P values when comparing individuals belonging to the highest quintile with individuals belonging to the lowest quintile, for the unweighted (OR = 1.36; 95% CI 1.07-1.73; P = 0.013) and weighted (OR = 1.39; 95% CI 1.08-1.80; P = 0.011) scores (Table 3).
We performed an additional analysis using all the individuals and the average scores rather than the absolute values (see Methods). We observed a strong association between long gdTL and increased risk of MM when analyzing the score as a categorical variable with a stronger association observed with the weighted score (OR = 1.69; 95% CI 1.36-2.11; P = 2.97 × 10 −6 for highest vs. lowest quintile) compared to the unweighted score (OR = 1.58; 95% CI 1.28-1.94; P = 1.93 × 10 −5 ) ( Table 3).
Considering that TERC-rs10936599 is a well-known risk locus for MM, we have also computed a score considering all the SNPs with the exception of this one, observing results that are in line with what we observed with the 11 SNP score (Supplementary Table 3).
We also estimated the correlation between MM risk and the genetic score using a Mendelian randomization

Association between individual SNPs and teloscore with MM OS
Survival analysis did not show any statistically significant associations if considering multiple testing. However, we found an association between TERC-rs10936599-T allele carriers and worse MM survival compared to C allele carriers (HR allelic = 1.20; 95% CI 1.00-1.43; P = 0.048) and an association between A allele homozygous of the PXK-rs6772228 SNP and worse survival compared to homozygous for the common T allele (HR homozygous = 4.51; 95% CI 1.06-19.17; P = 0.041) ( Table 4 and Supplementary Table 4).
The analyses between teloscore and OS of MM patients showed a trend between having longer gdTL and better MM survival, considering the quintile as a continuous variable for unweighted (HR = 0.93; 95%CI 0.85-0.99; P = 0.046) or weighted (HR = 0.93; 95% CI 0.86-0.99; The estimate measures the increase in risk associated with each increase of one quintile. P = 0.049) score, adjusted by DS staging system. Similar results were obtained in the analyses adjusted by ISS ( Table 5).
The results of the Mendelian randomization analysis performed with IVW regression and MR-Egger regression showed an inverse, non-significant association between MM survival and the genetic score built with the 11 SNPs (coeff = −0.44, 95% CI −1.18 to 0.29, P = 0.238), with the intercept of the MR-Egger regression non-significantly different from zero (coeff = −0.108, 95% CI −0.28 to 0.07, P = 0.227). Q′ heterogeneity (P = 0.391) did not show statistically significant heterogeneity between the polymorphisms, and Bowden I 2 GX statistic did not show weak instrument bias (92.8%). Scatter plots and forest plot of the causal estimate for risk of developing MM are shown in Supplementary Figs. 3 and 4, respectively.

In silico identification of functional effects
We tested the possible functional relevance of CTC1-rs3027234, ZNF208-rs8105767, OBFC1-rs9420907, TERC -rs10936599, and PXK-rs6772228, that were significant in at least one analysis and the variants in linkage disequilibrium (LD) with them using RegulomeDB and GTEx.
RegulomeDB showed a value of 1b for CTC1-rs3027234 that indicates the presence of an eQTL, of at least one transcription factor binding site and the region to be sensitive to the action of DNase. Analyzing rs3027234 in the GTEx portal, a correlation between the T allele and lower gene expression was observed in 20 tissues including whole blood, where it showed a strongly significant association (P = 1.2 × 10 −33 ). RegulomeDB assigned to ZNF208-rs8105767 a value of 1 f, which indicates a functional role similar to CTC1-rs3027234. GTEx reported six associations with different genes for rs8105767, and the strongest correlation was between the A allele and lower expression of the ZNF257 gene (P = 1.6 × 10 −25 ) in whole blood. RegulomeDB assigned to OBFC1-rs9420907 and TERC-rs10936599 a value of 3a, indicating the presence of at least one transcription factor binding site and the region to be sensitive to the action of DNase, and GTEx did not show any eQTL. RegulomeDB assigned to PXK-rs6772228 a value of 7, indicating the absence of a functional role. Despite that, rs6772228 is in LD with rs77480019 (r 2 = 0.92) and rs73077957 (r 2 = 0.84) for which RegulomeDB assigned a value of 3a, but GTEx did not show any eQTL in whole blood.

Discussion
Telomere length is considered a risk factor, or a risk marker, for a large number of diseases including several cancer types [18][19][20][21][22][23][24][25] . However, in spite of the number of studies, including several with considerable sample size, it is still unclear whether long or short telomeres are responsible for the increase in risk. In addition to true biological differences in various diseases, a possible explanation for this uncertainty may be explained by the technical heterogeneity due to sample handling and study design used across the studies. For MM only two relatively small studies have been attempted, one in a retrospective case-control setting 55 and the other in a prospective cohort 26 , both reporting an association between longer telomeres and increased risk of developing the disease.
Telomere length is at least in part genetically determined and there are several studies that have successfully used SNPs as genetic surrogates for telomere length to infer risk of developing several tumor types 23,[37][38][39][40][41][42][43][44][45] . This approach has never been attempted in MM. LTL has been tested in relation to cancer survival 27 , but never in MM and the only study conducted in MM is very small (n = 5) and used purified bone marrow cells and not whole blood 28 .
Considering that LTL measure could be influenced by several epidemiologic confounders or methodological issues, the main advantage of the strategy used in this study is that the exposure (the genotypes) is not influenced by study design (prospective vs. retrospective), by sample manipulation or by environmental variables since genotypes are invariable throughout life. The aim of this study was to test whether GWAS-identified SNPs involved in telomere length, combined in a score, were associated with MM risk and, for the first time, OS of MM patients. Considering that individually each of the SNPs explains a small proportion of LTL variability, we have combined them in a score and performed a MR analysis to better characterize the relation between gdTL and the disease risk and survival.
We found a novel association between carriers of the C allele OBFC1-rs9420907 SNP and increased MM risk (P = 7.78 × 10 −4 ). The C allele of the SNP is also associated with longer telomeres 34,38 . This gene product is associated to the shelterin complex and is also a subunit of alpha accessory factor (AAF) that is involved in the initiation of DNA replication. We observed an increase in risk associated with the C allele of the rs3027234 polymorphisms and with G allele of the rs8105767. The CTC1-rs3027234-C allele is associated with longer telomeres and, according to GTEx, with an increased expression of the gene in blood with a very high statistical significance (p = 1.2 × 10 −33 ). It is therefore plausible that the C allele contributes to an increase efficiency of the complex that in turn results in longer telomeres that increase the risk of developing MM. The G allele of rs8105767 is associated with longer telomeres 34 and, according to GTEX, it has an effect on the expression of various genes in whole blood, mainly ZNF257 (P = 1.6 × 10 −25 ).
The main finding of this manuscript is the observation that gdTL is associated with increased risk of MM. This association is consistent in all the models we tested and when using all the subjects it reached a very high statistical significance, especially comparing individuals with the longest telomeres with individuals with the shortest telomeres (quintile 5 vs. quintile 1; P = 2.97 × 10 −6 ) or when considering the quintile variable as continuous (P = 1.07 × 10 −8 ). This finding is in agreement with what was found in the two studies that reported longer telomeres measured by real-time PCR to be associated with increased risk of developing the disease 26,55 and it is also in line with a very large study performed on a substantial number of different cancer types, where the general trend for most cancers was an association between longer genetically determined telomeres and increased risk 23 . To analyze in depth the association between gdTL and risk of MM, the 11 SNPs selected were analyzed with a Mendelian randomization approach to test the presence of directional pleiotropy and its influence in the estimate of the causal effect between LTL and risk of developing MM. The results were in line with the finding obtained through teloscore analysis, supporting the evidence for a causal effect of longer gdTL with an increased MM risk. Moreover, the MR-Egger intercept was not statistically different from zero and the results of the heterogeneity tests were not statistically significant. These observations suggest the absence of directional pleiotropy and heterogeneity between the polymorphisms selected.
From a biologic point of view, longer telomeres (particularly as measured by the SNP teloscore, which represents the innate tendency to have longer or shorter telomeres) could represent a marker of a cell that has a higher dividing potential and therefore an increased tendency in acquiring new potential harmful mutations 23 .
A potential limitation is that most of the controls used in this study, as for almost all case-control studies, are either blood donors or hospitalized individuals with pathologies unrelated to MM. Considering that gdLTL has been associated with several human traits, this may represent a confounding factor of the analysis.
A point of novelty is represented by two borderline associations between the allele associated with shorter TL and a worse OS of MM patients, compared to the allele associated with longer TL, both for TERC-rs10936599 (P = 0.048) and PXK-6772228 (P = 0.041). Furthermore, the teloscore survival analyses showed that longer gdTL is correlated with better MM patients OS compared to patients with shorter gdTL. Despite this, the different direction obtained for IVW and MR-Egger analyses did not support a significant association of the teloscore with OS of MM patients, although the Cochran's Q statistic did not show evidence of heterogeneity between the polymorphisms. A possible explanation for Mendelian randomization results is that the increase of patient's survival could be due to pleiotropic effects of the genetic variants that act on the outcome via a confounder.
Xu and colleagues conducted a thorough meta-analysis to test the effect of LTL in cancer progression and survival and even though MM was not included, their findings are in agreement with ours i.e., shorter telomere associated with worse prognosis 27 .
We do not have data for myeloma-specific mortality. On the other hand, although in the last years many advances in therapy have been made, MM remains an incurable disease and therefore we presume that the vast majority of deaths are disease specific.
It is worth noting that gdTL as measured in our work, being based on germline polymorphisms that do not change throughout life, is unaffected not only by age, but also by any other biological mechanism induced by environmental exposures (e.g., smoking) or therapeutic exposures (e.g., anticancer treatments) that are known to affect telomere length.
In conclusion, we present here convincing evidence that longer gdTL is a risk marker for MM development.