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 cancer1. MM has usually a late onset with a mean age at diagnosis around 60 years1.

There are overwhelming evidences that genetic variability influences the risk of developing MM2,3,4,5,6,7,8. In the recent years the importance of genetics has also emerged in MM response to treatment and survival9,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 cells16. Telomerases are the most important enzymes involved in telomere replication and their dysfunction is considered a cancer hallmark17. Telomere length measured in leukocytes (LTL) has been used to investigate the susceptibility to different cancer types including several hematological malignancies18,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 study21,26. In a meta-analysis conducted by Xu and colleagues the role of LTL in relation to survival in several cancer types has been investigated27. 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 cells28.

Telomere length is highly correlated across tissues29,30, it is reasonably stable during time, as assessed in a longitudinal study18, 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.7032,33. Recently genome-wide association studies (GWAS) have identified 11 SNPs23,34,35,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 types23,37,38,39,40,41,42,43,44,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.

Materials and methods

Study population

This study was conducted within the IMMeNSE consortium, that has been extensively described elsewhere46. 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) criteria47. In addition, clinical data such as disease stage according to Durie-Salmon (DS) staging system48 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 IMMEnSE study protocol was approved by the Ethics Committee of the Medical Faculty of the University of Heidelberg (reference number: S-004/2020). Following the guidelines of the Declaration of Helsinki, written informed consent was obtained from each participant.

Table 1 Summary of the study population.
Full size table

SNP selection

To build the teloscore we used 11 single nucleotide polymorphisms (SNPs) (ZNF676-rs412658, TERT-rs2736100, CTC1-rs3027234, DHX35-rs6028466, PXK-rs6772228, NAF1-rs7675998, ZNF208-rs8105767, OBFC1-rs9420907, ACYP2-rs11125529, TERC-rs10936599, and ZBTB46-rs755017), identified through GWAS, that affect telomere length23,34,35,36. Characteristics of the SNPs under investigation are summarized in Supplementary Table 1.

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 elsewhere37. 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 estimates50,51. The Egger (MR-Egger) estimation method was used to test for possible pleiotropic effects of the genetic instruments50,51. Finally to assess evidence of heterogeneity between the polymorphisms, the Cochran Q′ heterogeneity and Bowden I2GX statistics were computed to obtain an estimate of instrument strength for the MR-Egger analysis52. Additional information has been given elsewhere50,51. Analyses were carried out with the Mendelian randomization package in R V0.5.0 (30 Sept. 2020).

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. The SNPs under investigation are summarized in Supplementary Table 1. All the selected SNPs were in Hardy–Weinberg Equilibrium (HWE) in the control population. The average call rate of the SNPs used for statistical analysis was 96.3%. Concordance with duplicated samples was higher than 99%.

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 (ORhomozygous = 0.55; 95% CI 0.39–0.76; P = 2.87 × 10−4), while OBFC1-rs9420907-C was associated with increased MM risk (ORheterozygous = 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.

Table 2 Association between single SNPs and MM risk.
Full size table

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).

Table 3 Association between teloscore and MM risk.
Full size table

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 approach as previously described. The IVW regression showed a highly significant association between the 11 SNPs used as genetic instruments of exposure and risk of developing MM (coeff = 1.175, 95% CI 0.61–1.75, P = 5.38 × 10−5), whereas the MR-Egger regression showed a statistically non-significant association (coeff = 0.532, 95% CI −1.16 to 2.22, P = 0.538), although in line with IVW regression, with the intercept slightly, but not significantly, different from zero (intercept = 0.06, 95% CI −0.09 to 0.20, P = 0.428). The Q′ heterogeneity statistic did not show statistically significant heterogeneity between the polymorphisms (P = 0.111), and Bowden I2GX statistic showed that the selected genetic variants have high power to detect directional pleiotropy and causal effect for MR-Egger (I2 = 93.1%). Scatter plots and forest plot of the causal estimate for risk of developing MM are shown in Supplementary Figs. 1 and 2, respectively.

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 (HRallelic = 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 (HRhomozygous = 4.51; 95% CI 1.06–19.17; P = 0.041) (Table 4 and Supplementary Table 4).

Table 4 Association between individual SNPs and MM OS, adjusted by stage (Durie-Salmon).
Full size table

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; P = 0.049) score, adjusted by DS staging system. Similar results were obtained in the analyses adjusted by ISS (Table 5).

Table 5 Association between teloscore and MM OS.
Full size table

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 I2GX 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 (r2 = 0.92) and rs73077957 (r2 = 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 types18,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 setting55 and the other in a prospective cohort26, 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 types23,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 survival27, 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 blood28.

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 telomeres34,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 telomeres34 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 disease26,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 risk23. 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 mutations23.

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 prognosis27.

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.