- Open Access
A variability in response of osteoclasts to zoledronic acid is mediated by smoking-associated modification in the DNA methylome
Clinical Epigenetics volume 15, Article number: 42 (2023)
Clinical trials have shown zoledronic acid as a potent bisphosphonate in preventing bone loss, but with varying potency between patients. Human osteoclasts ex vivo reportedly displayed a variable sensitivity to zoledronic acid > 200-fold, determined by the half-maximal inhibitory concentration (IC50), with cigarette smoking as one of the reported contributors to this variation. To reveal the molecular basis of the smoking-mediated variation on treatment sensitivity, we performed a DNA methylome profiling on whole blood cells from 34 healthy female blood donors. Multiple regression models were fitted to associate DNA methylation with ex vivo determined IC50 values, smoking, and their interaction adjusting for age and cell compositions.
We identified 59 CpGs displaying genome-wide significance (p < 1e−08) with a false discovery rate (FDR) < 0.05 for the smoking-dependent association with IC50. Among them, 3 CpGs have p < 1e−08 and FDR < 2e−03. By comparing with genome-wide association studies, 15 significant CpGs were locally enriched (within < 50,000 bp) by SNPs associated with bone and body size measures. Furthermore, through a replication analysis using data from a published multi-omics association study on bone mineral density (BMD), we could validate that 29 out of the 59 CpGs were in close vicinity of genomic sites significantly associated with BMD. Gene Ontology (GO) analysis on genes linked to the 59 CpGs displaying smoking-dependent association with IC50, detected 18 significant GO terms including cation:cation antiporter activity, extracellular matrix conferring tensile strength, ligand–gated ion channel activity, etc.
Our results suggest that smoking mediates individual sensitivity to zoledronic acid treatment through epigenetic regulation. Our novel findings could have important clinical implications since DNA methylation analysis may enable personalized zoledronic acid treatment.
Bisphosphonates are used to control and reduce bone loss in a variety of bone pathological conditions, such as cancer-induced bone disease and osteoporosis [1, 2].
Bisphosphonates belong to a class of drugs that are targeted specifically to mineralized bone due to its unique ability, binding to hydroxyapatite. Upon bone resorption, the bone resorbing osteoclast takes up the drug. At this stage, the side-chain of the bisphosphonates will trigger an inhibition of the bone resorptive activity. The N-containing bisphosphonates are reported to target the mevalonate pathway through inhibition of the farnesyl diphosphate synthase (FDPS) (for more details please refer to [1, 2]). Zoledronic acid is a N-containing bisphosphonate and is one of the most potent bisphosphonates to prevent bone loss, as demonstrated in several clinical trials [3, 4]. Despite this, it is also known that it is not equally potent on all patients [5,6,7,8,9]. In the case of cancer patients with bone metastasis around 50% of patients treated with zoledronic acid develop new skeletal-related events despite treatment for 1 year, while with placebo treatment this rate is 70% [5, 6]. Thus, zoledronic acid treatment is not potent on all cancer patients with bone disease. Also for the treatment of osteoporosis, variations in potency of zoledronic acid treatment between individuals are observed [3, 9]. Why may there be such a variation in potency from one person to another? Without a doubt, an answer to this will be multi-factorial, but one priority would be to understand if it is influenced by environmental factors including life-style.
It is well known that the dynamic relationship between bone resorption and formation, which is needed to maintain bone mass and bone health, can be influenced by multiple factors. These include biological factors like age, gender, and menopausal status , but also life-style factors like alcohol overuse [11, 12] and tobacco smoking [12,13,14,15]. However, it is unclear whether these factors can also affect the potency of drug treatment when using, e.g. zoledronic acid, it seems that insufficient specific knowledge is available regarding this issue . Recent studies suggest that an inadequate response of osteoporosis patients to primarily bisphosphonates correlates with the presence of resistant osteoclasts . This could suggest that an incomplete response of patients to zoledronic acid could be due to a variation in sensitivity of osteoclasts.
We have recently investigated the sensitivity of osteoclasts in vitro to zoledronic acid. These osteoclasts were generated from peripheral blood CD14+ monocytes obtained from 46 healthy women. We identified the half-maximal inhibitory concentration (IC50) of zoledronic acid based on total eroded bone surface for each of the osteoclast donors and found a surprising > 200-fold inter-individual variation in sensitivity. Our analyses revealed that smoking was a significant contributor to this variation . Given that smoking is well known to have strong influences on the epigenetic regulation of genes, in particular DNA methylation [19, 20], it is reasonable to speculate that smoking may affect the sensitivity to zoledronic acid through epigenetic regulation. Previous epigenome-wide association studies (EWASs) have demonstrated that cigarette smoking reduces DNA methylation levels at multiple genomic loci in blood cells [21, 22]. Both prenatal and current cigarette smoke exposures have been associated with reduced DNA methylation in genes involved in chemical detoxification such as CYP1A1 and even hypomethylation at genome level [23, 24]. It is well realized that the environment modulates genetic effects , although human genomic studies, e.g. genome-wide association studies (GWASs), and expression quantitative trait loci (eQTL) studies rarely test for genetic interactions with environmental exposures.
Epigenetics refer to the meta-level regulation of gene expression caused by mechanisms other than changes in the DNA sequence. Under the constant influence of external factors, epigenetic mechanisms regulate which genes are turned on and off to adapt their expression to a change in environment. As such, epigenetics serves as the bridge between nature, our genome, and nurture, our environment . Different molecular mechanisms are involved in epigenetics, including DNA methylation, histone modifications, and non-coding RNAs (ncRNAs) [28, 29]. In this regard, epigenetic regulation through DNA methylation involves the transfer of a methyl group onto the C5 position of the cytosine to modify the function of the DNA, e.g. to result in gene silencing . In recent studies, we have reported that monocytes, precursors of the human osteoclasts, could have been epigenetically programmed through DNA methylation to reflect age and menopausal status of women, resulting in more aggressive osteoclasts [30, 31]. The significant associations of environmental, including behaviour, factors with bone metabolism call for epigenetic studies to elucidate the molecular mechanisms underlying the interplay between the exogenous factors, nurture, and the genome, nature.
Therefore, we have investigated the smoking-induced epigenetic changes in the genome of peripheral blood mononuclear cells and the possible influence of these epigenetic changes on the sensitivity of resulting osteoclasts to zoledronic acid. We performed an epigenome-wide multifactorial association study on a cohort of pre- and post-menopausal women to investigate the differential DNA methylation regulation patterns associated with the sensitivity to zoledronic acid while taking into consideration the number of cigarettes they have smoked throughout life. Our analyses reveal a strong epigenome-wide association of the DNA methylation patterns at 59 CpG sites with the smoking-dependent sensitivity of osteoclastic bone resorption to zoledronic acid treatment.
Characteristics of blood donors
The study group used for the current analyses consists of 34 blood donors out of the 46 used in our previous publication by Møller et al. . The donor characteristics of the original study group can be seen in the original publication , while the characteristics of the current sub-group can be seen in Table 1. There are no statistically significant differences between the original samples and the sub-group of samples used for this study with respect to any of the demographic characteristics (mean age 52 vs 53 years—t-test p = 0.522; median weight 73 vs 72.5 kg—Mann–Whitney test p = 0.956; mean height 169.5 vs 169.7 cm—t-test p = 0.902). In our previous study, we only had information on current smoking status, but after the telephone interview, we can conclude that 12 had never smoked, 16 were past smokers, and 6 were current smokers (Table 1). The median number of cigarettes smoked through a lifetime, including non-smokers, was 6096, but with a large range from 0 to 268,800 (Table 1). The IC50s of zoledronic acid on osteoclastic bone resorption in vitro spanned from 0.061 to 9.49 µM reflecting a 156-fold difference from min to max and with a median of 0.28 µM (Table 1). A comparison of the IC50 values between the original and present study group showed that they did not differ, median 0.26 µM and 0.28 µM, respectively (Mann–Whitney test p = 0.8332). Hence, the current sub-group of 34 does not differ from the original study group of 46.
EWAS on zoledronic acid IC50 of osteoclast cultures and smoking
The EWAS with 865,857 CpG sites was performed using the regression model with DNA methylation as a dependent variable, IC50 and smoking as main effect explanatory variables together with their interaction for assessing the smoking-dependent association between IC50 and DNA methylation. The analysis was corrected for age and cell composition effects. For the smoking-dependent association of IC50 with DNA methylation (i.e., the interaction effect), we identified 59 CpGs displaying genome-wide significance with FDR < 0.05 (corresponding p-value < 4.27e−06) (Table 2, Additional file 3: Table S1). A Manhattan plot for the p-values and chromosomal location of all the CpGs can be seen in Fig. 1. The QQ plot in Fig. 2 also exhibits a clear pattern of deviation from being random starting roughly from CpGs with p < 1e−03. Among the top significant sites, 3 CpGs have p < 1–08 and FDR = 2e−03 (cg00227784, cg14355428, cg22010000) (Table 2). The QQ plot in Fig. 2 shows that the top significant sites for smoking-dependent association with IC50 deviate remarkably from the diagonal line of random association. In the volcano plot in Fig. 3, all CpGs reaching genome-wide significance are marked with red color and the CpGs showing suggestive significance are marked as blue. Among the red colored top CpGs, we could see more CpGs with decreased methylation as compared to CpGs with increased methylation.
We identified 1 CpG of genome-wide significance for association with the main effect of IC50 (cg00133289, FDR = 0.013, p = 1.91e−08) (Additional file 3: Table S1). This CpG relates to the tripartite motif-containing 22 gene (TRIM22), coding for a nuclear protein that functions as a nuclear E3 ubiquitin ligase. Functional annotation of differentially methylated CpGs was performed by over-representation analysis of gene ontology (GO) terms using clusterprofiler . Over-representation analysis was performed on 3934 unique genes linked to the 21,760 CpGs associated with the interaction effect between IC50 and smoking with p < 0.05. We found 18 GO clusters significantly over-represented by the 3934 genes with FDR < 0.05 (Table 3, Fig. 4). Among the top significant sites are protein serine kinase activity, frizzled binding, minus-end-directed microtubule motor activity, gated channel activity, cation:cation antiporter activity, chloride channel regulator activity, protein serine/threonine kinase activity, signaling adaptor activity, etc.
Enrichment analyses for the 59 CpG regions showed that at distances ranging from 1,000 to 50,000 bp to the UK Biobank GWAS signals, there was a significant enrichment of 2 genes related to GWAS on bone traits. In addition, 14 genes related to GWAS on body size measures when compared to unrelated GWAS on mental disorders (Table 4; Additional file 1: Figure S1 and Additional file 2: Figure S2).
As an independent validation, we further analyzed the top 59 genome-wide significant CpGs displaying genome-wide significance (p < 1e−08) with a FDR < 0.05 for the smoking-dependent association with IC50, in a dataset  with multi-omics analyses of BMD in 119 Caucasian female subjects. As indicated in Table 2, among our 59 CpG sites with the smoking-dependent sensitivity to zoledronic acid treatment, 29 were in close vicinity (within 5 kb up- and downstream) of genomic sites significantly associated with BMD (while 30 CpGs had no associations with q-value < 0.05 in the dataset ). Remarkably, cg04932413 on chr. 14, cg00013660 on chr. 9, cg22581270 in CASP16P, and cg14485214 in SNED1 were highly significantly associated with BMD (p < 7.0e−07, FDR < 4.7e−05) (Additional file 4: Table S2). For the main effect of 17 CpGs on IC50 with FDR < 0.1 (Additional file 3: Table S1), 9 CpGs were validated with significant association with BMD (Additional file 5: Table S3).
We have performed EWAS to identify DNA methylation patterns associated with the sensitivity (IC50) of osteoclasts from different individuals to zoledronic acid as well as whether these associations are smoking-dependent. We found 59 CpGs reaching genome-wide significance for the smoking-mediated association with IC50 together with three for the main effects of IC50.
Zoledronic acid, a bisphosphonate, is commonly used to treat osteoporosis and bone metastases. Bisphosphonates target the osteoclast in order to prevent bone resorption [1, 2] and are in general very potent to prevent bone loss, as shown in a multitude of clinical studies [3,4,5,6]. However, while the majority of patients have clear benefits of bisphosphonate treatment and seem to be responsive to this treatment, it is not so for all patients. Cairoli et al.  conducted a study where 97 patients, who started treatment for osteoporosis with bisphosphonates, alendronate or risedronate, were followed for 3 years. The authors found that 25 out of the 97 patients matched the criteria of treatment failure  by meeting at least one of the following criteria: (i) two or more incident fragility fractures and (ii) a decrease in BMD greater than the least significant change (LSC) . In their study, Cairoli et al.  found two variables that independently correlated with treatment failure, namely smoking and elevated plasma levels of alkaline phosphatase (ALP) at baseline. Current smoking alone displayed an odds ratio of 3.2, while smoking combined with ALP-levels reached 8.03. This strongly suggests that current smoking increases the risk of treatment failure when using risedronate or alendronate to treat osteoporosis. Only few other studies have addressed the issue and they did not observe any consequence of smoking on the potency of bisphosphonates [16, 36, 37]. In general, these studies primarily focused on the current smoking status or combined it with past smoking. This rough definition of smoking will likely oversimplify something that is very individual, namely the smoking history. Therefore, a strength of our study (though small and conducted ex vivo) is the detailed information on the individual number of cigarettes smoked through life.
In our recent study by Møller et al. , current smoking was found to significantly correlate with the sensitivity (IC50) of human osteoclasts prepared from different donors. This phenomenon could suggest smoking as a major contributor to the more than 200-fold variation in sensitivity to zoledronic acid observed between the osteoclasts from different donors . For the current study, new telephone interviews with the participants allowed us to obtain more details on the number of cigarettes smoked through life. This gave us the opportunity to better investigate the potential influence of smoking on the sensitivity of osteoclasts to zoledronic acid.
This updated information allowed us to do EWAS analyses where we investigated if there is a DNA methylation profile of blood cells that may predict the sensitivity (IC50) of osteoclasts to zoledronic acid in the context of cigarettes smoked throughout life. Our analyses indeed identified 59 unique CpGs where the DNA methylation level reflected the number of cigarettes smoked interacting with the IC50 of zoledronic acid. These could be linked to 37 genes. It is important to stress that these 59 CpGs are identified in the DNA based on a blood sample. Since the sensitivity to zoledronic acid is determined in the donor-derived osteoclasts after 10 days of differentiation, it means that the methylation status of the identified CpGs is persistent and may have direct consequences for the sensitivity of the resulting osteoclasts.
In order to identify, which cellular processes are represented by the 21,760 CpGs nominally-significant (p < 0.05) for the smoking-dependent association with IC50, we analyzed the 3934 genes linked to these CpGs for over-representation of gene ontology (GO) terms. The analysis identified 18 enriched GO biological processes (FDR < 0.05). It is striking that just a few types of cellular processes are enriched: actin and microtubule-related functions, membrane related transporters and channels, etc. These processes match well-known targets or consequences of N-containing bisphosphonates [38,39,40,41,42,43,44,45,46] suggesting why the functions and pathways listed in Fig. 4 and Table 3 may affect osteoclasts’ sensitivity to zoledronic acid.
To investigate further the possible function of the 59 CpGs identified in the interaction analysis, we performed an enrichment analysis with existing GWAS databases on bone mineral density and body size measures compared to the unrelated phenotypes—mental disorders. We illustrated the impact of distance from a CpG on the enrichment for GWAS-significant loci associated with relevant traits. By doing so, we identified 15 genes that were significantly enriched within a ≤ 50,000 bp distance between the CpG and SNP for either bone mineral density or body size. These genes are likely to relate to the sensitivity of osteoclasts to zoledronic acid.
In our replication study based on the dataset of Qiu et al. , we found that 29 out of the 59 significant CpGs were within 5 k bp of CpGs associated with BMD. Of these 29 CpGs, 8 were hyper- and 21 were hypomethylated, suggesting that most of their related genes will be expressed at higher levels and that this would render osteoclasts more resistant to zoledronic acid. Interestingly, nine of the CpGs found to be significantly associated with BMD (Table 2, Additional file 4: Table S2), were also identified in the GWAS enrichment analysis (shown in Table 4). This therefore suggests a more direct link between these nine genes, their predicted expression level based on DNA methylation changes due to smoking, and BMD as well as sensitivity to zoledronic acid. It is generally known that reduced DNA methylation in the promoter region activates gene expression, while in the intragenic region it may rather reduce expression or even modulate alternative splicing [47, 48]. All nine CpGs were found to be demethylated by smoking, suggesting an increased (if located in promoter) or possibly decreased/aberrant (if located in gene body) expression of related genes. The nine CpGs and their connected genes are: cg15073625 (KLHL29, body), cg17641218 (VAMP8, promoter), cg06342954 (AGAP1, body), cg14485214 (SNED1, body), cg06501109 (DDR1, promoter), cg14355428 (CYP51A1, body), cg08086799 (TTC9C, promoter), cg01279902 (XRCC3, body), and cg04662961 (FRMD5, body).
The aim of our study was to reveal the molecular basis for how DNA methylation changes in peripheral blood mononuclear cells, triggered by smoking, may cause osteoclasts to become less sensitive to zoledronic acid treatment. Based on all of our presented data, the nine listed CpGs may serve as candidate markers for the underlying molecular basis. However, not much is known about these in the context of sensitivity to zoledronic acid or bisphosphonates. Yet, the following may be of special interest: (1) Vesicle-associated membrane protein 8 (VAMP8), belongs to the SNARE family facilitating membrane fusion between late endosomes as well as in autophagy and is closely linked to VAMP7 [49, 50]. Endosomes play a key role in enabling zoledronic acid to be released into the cytosol of the osteoclast and to act on its target . (2) CYP51A1 has, to our knowledge, not been reported to have a function specifically in bone cells or osteoclasts. CYP51A1 is known to be involved in the biosynthesis of cholesterol and acts downstream of FDPS in the mevalonate pathway . However, since we are addressing drug-related effects on the mevalonate pathway, which is not specific to osteoclasts, other non-bone related pathways may also come into play. CYP51A1 is therefore active in the same pathway that is targeted by zoledronic acid and may therefore explain why a CpG, regulating this particular gene, is the second most significant in our EWAS analysis (Table 2). Furthermore, patients with mutations in CYP51A1, amongst other things, display bone defects and reduced cholesterol levels . Thus, the significant validations in Additional file 4: Table S2 and Additional file 5: Table S3 provide further evidence that our identified top associated sites are potentially under regulation by meQTLs, with DNA methylation as an important epigenetic mechanism. This may mediate the interaction between environmental exposure (here smoking) and the genome regulation.
We acknowledge that our study has limitations. A major limitation is the number of samples. More samples may have resulted in additional genome wide significant CpGs, in particular for CpGs related to IC50 alone. This would have allowed a more general view on what determines the sensitivity of osteoclasts to zoledronic acid. An additional limitation is that a comparative analysis should be performed for the identified CpGs in terms of gene expression evaluation, something that has recently been initiated. Finally, our study was done using osteoclasts ex vivo and this of course limits the level of interpretation and impact. However, our study also has multiple strengths. We have foremost used primary human cells as opposed to cells from animals or cell lines. The use of osteoclasts ex vivo as the experimental model allows us to identify specific genes that directly or indirectly play a key role in the sensitivity of the target cells to zoledronic acid. This would not have been possible if only using clinical data from patients undergoing treatment. The combination of detailed demographic information, detailed experimental analysis, and epigenetic analysis has allowed us to obtain strong EWAS data showing that smoking epigenetically regulates genes that are linked to the sensitivity of osteoclasts to zoledronic acid.
We have identified a mechanism that may explain why past and present smoking negatively affects the clinical potency of drugs such as zoledronic acid . Therefore, this study goes beyond merely identifying changes of DNA methylation due to smoking, since we specifically identify those CpGs that are affected by smoking and affect the downstream drug sensitivity of osteoclasts. This has not been shown before. It may give us the possibility to employ a DNA methylation profile of peripheral blood mononuclear cells as a clinical tool to determine who may have optimal or limited benefit of treatment with zoledronic acid. This could enable a more personalized approach to treatment of patients with respect to choice of drug, dosing, and duration of treatment. Of course, this would demand clinical testing to validate such an approach; such investigations are currently ongoing. However, a follow-up to this study by experimental settings should also be conducted to inform the functionality of our candidate meQTLs and further confirm causal genes.
Material and methods
The study samples
The study sample consists of 34 female blood donors aged between 40 and 66 (Table 1). They were recruited for the voluntary participation as regular blood donors at the Bloodbank of Lillebaelt Hospital, Denmark, based on the ethical approval on the 11th of May 2015 by the Danish Regional Scientific Ethics Committee (S-20150059) and signed informed consent. A questionnaire was used to collect their personal characteristics, medical history, and lifestyle, while a telephone interview was used to obtain detailed information on their past and present history of cigarette smoking. Through this interview, the following information was collected: (1) present smoker, (2) past smoker, (3) smoked for how many years, (4) stopped smoking when, (5) cigarettes smoked per day. Based on this information we could estimate the total number of cigarettes smoked through their lifetime. The 34 participants reflect a sub-fraction of the participants already used for other publications [18, 30, 31]. These 34 reflect the number of participants that we were able to re-contact for the interview and with respect to demographics and experimental out-come they are fully representative of the original 46.
Blood sample collection
From each donor, the buffy coat from approximately 500 ml of blood was used to isolate CD14+ monocytes for generation of osteoclasts. Furthermore, two 4 ml blood samples were collected when the donors were in a fasting state early in the morning 14 days later than their 500 ml blood donation. One sample was used for extraction of DNA from whole blood. The 4 ml blood samples for collecting DNA were stored at -80 °C.
Generation of osteoclasts
Osteoclasts were generated according to a standard procedure as described in . In brief, CD14+ monocytes were isolated, exposed to 25 ng/ml M-CSF (R&D System, Abingdon, UK) for 2 days followed by exposure to 25 ng/ml of both M-CSF and RANKL (R&D System) for another 7 days with renewal of media twice. Cells were cultured in αMEM (Invitrogen, Carlsbad, CA, USA), 10% FBS (Sigma-Aldrich, St. Louis, MO, USA), and 5% CO2 at 37 °C in a humidified atmosphere.
Determination of osteoclast sensitivity to zoledronic acid and their IC50
In order to determine the sensitivity of osteoclasts of individual donors in vitro the differentiated mature osteoclasts were lifted by Accutase (Biowest BW, Nuaillé, France) treatment and reseeded onto bovine cortical bone slices (BoneSlices.com, Jelling, Denmark), pre-coated with different concentrations of zoledronic acid. Data for the present study was obtained from our previously published study and a detailed description of the procedure can therefore be found in .
DNA methylation analysis
DNA methylation analysis was performed according to the Infinium HD Methylation Assay Reference Guide (document # 15,019,519 v07; https://support.illumina.com/downloads/infinium-hd-methylation-reference-guide-15019519.html). Genomic DNA (0.2–1.0 μg) was bisulfite converted using an EZ DNA Methylation-Direct Kit (Zymo Research, Irvine, CA, USA). DNA samples were bisulfite converted by incubation with the CT conversion reagent for 8 min at 98 °C, 3.5 h at 64 °C, followed by 18 h at 4 °C in a thermocycler. The treated DNA was added to a Zymo-Spin IC Column, desulfonated using M-desulfonation buffer, and then eluted from the column in 12 μl of M-elution buffer. Methylation profiling of the bisulfite-treated DNA was performed using Illumina Infinium MethylationEPIC BeadChip (Illumina, San Diego, CA, USA) according to standard protocol. In brief, 4 μl of bisulfite-treated DNA was denatured, neutralized and amplified with an overnight whole-genome amplification reaction. The amplified DNA was then enzymatically fragmented, precipitated and re-suspended in hybridization buffer before being dispensed onto the MethylationEPIC BeadChips for hybridization. After hybridization, the BeadChips were processed through a primer-extension protocol and subsequently stained. Finally, the BeadChips were coated and imaged using Illumina’s HiScan System.
DNA methylation data pre-processing
The raw methylation data cover 865,857 CpG sites across the genome. The R package minfi  was used for data pre-processing including quality control (QC) and normalization. For each CpG site of a sample, a detection p-value was first calculated by comparing the total signal for each probe to the background signal level estimated from the negative control probes. Very small p-values are indicative of a reliable signal whilst large p-values generally indicate a poor quality signal. We filtered out 2221 CpGs with detection p-value > 0.01 in more than 5% of the overall samples (i.e. 2 samples). After QC, a total of 863,686 CpG sites remained. We further removed all CpGs on the Y-chromosome (535 CpGs) and CpGs physically overlapping with SNPs (181,195 CpGs) leaving 683,408 CpGs for subsequent analysis. We kept X-linked CpGs because all our samples are females. QC at sample level was done by plotting the log median intensity in the methylated (M) against that of the unmethylated (U) channels using getQC and plotQC functions in minfi, with no bad quality sample found. Data normalization was performed by the functional normalization  implemented in the minfi R-package. At each CpG site, DNA methylation level was summarized by calculating a methylation “beta”-value defined by the Illumina’s formula as β = M/(M + U + 100). Before statistical analysis, the β-values were converted to methylation M-values for better statistical properties by logit transformation with M = log2(β/(1-β)) .
Controlling cell-type composition
Since the target tissue is whole blood comprising multiple cell types, cellular heterogeneity among samples can be an important confounding factor in epigenetic association analysis due to cell specificity of DNA methylation. To control for cell-type composition effect, we introduced ReFACTor, a reference-free adjustment for cell-type composition based on principal component analysis (PCA) . The algorithm calculates components that are correlated with the cell-type composition of the samples by applying an unsupervised feature selection step followed by PCA. Instead of estimating absolute cell count values, ReFACTor calculates the linear transformations of the cell-type composition as PCA components.
We validated the top associated CpGs displaying genome-wide significance (p < 1e−08) with a FDR < 0.05 for the smoking-dependent association with IC50, using an independent multi-omics study for BMD , which contains DNA methylome data from peripheral blood monocytes in 119 Caucasian female subjects. The selected CpG sites identified in this study, which are located in close vicinity (within 5 kb up- and downstream) with CpGs significantly associated with BMD (FDR q-value < 0.05) in the validation study were considered as validated.
We applied the hypergeometric test for over-representation analysis (ORA) to assess if the overlap of identified markers with those from a functional cluster or category (e.g., biological pathway) is significantly different from being random by calculating a probability from the hypergeometric distribution. ORA has also been implemented in a R package for biological pathway analysis, clusterProfiler , to test if members of one biological pathway are over-represented in a list of identified genes.
Enrichment of GWAS associations in vicinity of the CpGs
We tested the proximity between identified significant CpGs and GWAS-reported SNPs within a specified distance (in base pairs). As stated above, CpGs physically overlapping with SNPs (181,195 CpGs) were filtered out. We downloaded UK Biobank GWAS summary statistics made available by the Neale lab (initially made public on August 1, 2018; http://www.nealelab.is/uk-biobank/). Specifically, we downloaded the sex-combined GWAS summary statistic files listed in the https://docs.google.com/spreadsheets/d/1kvPoupSzsSFBNSztMzl04xMoSC3Kcx3CrjVf4yBmESU/edit#gid=227859291. We kept only traits marked as ‘high confidence’, with estimated heritability > 0.01, and h2_z (Z-scores for test of h2 > 0) ≥ 7. After we filtered out SNPs with p-value < 5e−6 there was a total of 393 GWAS datasets.
We further selected the GWAS for traits deemed “relevant” and less relevant for bone physiology. The following grouping was tested: Phen1: BMD and related traits (e.g. arthritis, osteoporosis treatments); Phen2: Body size (weight, height, BMI, length of body parts etc.); Phen3: Smoking and related traits; Phen4: Mental disorders as a “negative control”. The number and list of phenotypes in each ‘Phen#’ group is provided in (Additional file 6: Table S4).
For each CpG and for each Phen#, success ratio was calculated as a proportion of #successful SNPs for that CpG to #all SNPs associated with that phenotype group, Phen#. “Successful” means that the distance between that SNP and particular CpG is less than the stated maximal distance. Series of incremental distances from 1000 to 50,000 bp were tested, and then compared across the phenotype groups. The one-tailed Mann–Whitney U test was used to compare differences in success ratios between Phen1 and Phen 2 and “negative control” Phen4 (Additional file 7: Table S5).
We applied the linear regression models to detect the association of DNA methylation with IC50 and smoking as main effects and their interaction effect, adjusting for age and cell composition:
Smoking was defined as the total number of cigarettes smoked life-long; age was defined as age at blood sampling. Considering the limited sample size of the study, we selected the top 3 principal components (PCs) to add as covariates in the regression model to control for cell-type heterogeneity.
To control for multiple testing, we calculated the false discovery rate (FDR) following Benjamini et al.  using the R function p.adjust. We define p < 1e−05 as suggestive significant and FDR < 0.05 as genome-wide significance.
Availability of data and materials
The datasets used and/or analyzed during the current study are available from the corresponding author on reasonable request.
Russell RG. Bisphosphonates: the first 40 years. Bone. 2011;49(1):2–19.
Cremers S, Drake MT, Ebetino FH, Bilezikian JP, Russell RGG. Pharmacology of bisphosphonates. Br J Clin Pharmacol. 2019;85(6):1052–62.
Langdahl BL. Overview of treatment approaches to osteoporosis. Br J Pharmacol. 2021;178(9):1891–906.
Ibrahim A, Scher N, Williams G, Sridhara R, Li N, Chen G, et al. Approval summary for zoledronic acid for treatment of multiple myeloma and cancer bone metastases. Clin Cancer Res. 2003;9(7):2394–9.
Lipton A, Fizazi K, Stopeck AT, Henry DH, Brown JE, Yardley DA, et al. Superiority of denosumab to zoledronic acid for prevention of skeletal-related events: a combined analysis of 3 pivotal, randomised, phase 3 trials. Eur J Cancer. 2012;48(16):3082–92.
Rosen LS, Gordon D, Tchekmedyian NS, Yanagihara R, Hirsh V, Krzakowski M, et al. Long-term efficacy and safety of zoledronic acid in the treatment of skeletal metastases in patients with nonsmall cell lung carcinoma and other solid tumors: a randomized, phase III, double-blind, placebo-controlled trial. Cancer. 2004;100(12):2613–21.
Lipton A, Smith MR, Fizazi K, Stopeck AT, Henry D, Brown JE, et al. Changes in bone turnover marker levels and clinical outcomes in patients with advanced cancer and bone metastases treated with bone antiresorptive agents. Clin Cancer Res. 2016;22(23):5713–21.
MacLean C, Newberry S, Maglione M, McMahon M, Ranganath V, Suttorp M, et al. Systematic review: comparative effectiveness of treatments to prevent fractures in men and women with low bone density or osteoporosis. Ann Intern Med. 2008;148(3):197–213.
Diez-Perez A, Adachi JD, Agnusdei D, Bilezikian JP, Compston JE, Cummings SR, et al. Treatment failure in osteoporosis. Osteoporos Int. 2012;23(12):2769–74.
Wheater G, Elshahaly M, Tuck SP, Datta HK, van Laar JM. The clinical utility of bone marker measurements in osteoporosis. J Transl Med. 2013;11:201.
Sampson HW. Alcohol’s harmful effects on bone. Alcohol Health Res World. 1998;22(3):190–4.
Prieto-Alhambra D, Turkiewicz A, Reyes C, Timpka S, Rosengren B, Englund M. Smoking and alcohol intake but not muscle strength in young men increase fracture risk at middle age: a cohort study linked to the Swedish national patient registry. J Bone Miner Res. 2020;35(3):498–504.
Al-Bashaireh AM, Haddad LG, Weaver M, Chengguo X, Kelly DL, Yoon S. The effect of tobacco smoking on bone mass: an overview of pathophysiologic mechanisms. J Osteoporos. 2018;2018:1206235.
Du Y, Li P, Wen Y, Liang X, Liu L, Cheng B, et al. Evaluating the correlations between osteoporosis and lifestyle-related factors using transcriptome-wide association study. Calcif Tissue Int. 2020;106(3):256–63.
Li H, Wallin M, Barregard L, Sallsten G, Lundh T, Ohlsson C, et al. Smoking-induced risk of osteoporosis is partly mediated by cadmium from tobacco smoke: the MrOS Sweden study. J Bone Miner Res. 2020;35(8):1424–9.
Eastell R, Black DM, Boonen S, Adami S, Felsenberg D, Lippuner K, et al. Effect of once-yearly zoledronic acid five milligrams on fracture risk and change in femoral neck bone mineral density. J Clin Endocrinol Metab. 2009;94(9):3215–25.
Leger B, Fardellone P, Cormier C, Ostertag A, Funck-Brentano T, Fabre S, et al. Inadequate response to treatment reveals persistent osteoclast bone resorption in osteoporotic patients. Bone. 2021;153:116167.
Moller AMJ, Delaisse JM, Olesen JB, Bechmann T, Madsen JS, Soe K. Zoledronic acid is not equally potent on osteoclasts generated from different individuals. JBMR Plus. 2020;4(11):e10412.
Gao X, Jia M, Zhang Y, Breitling LP, Brenner H. DNA methylation changes of whole blood cells in response to active smoking exposure in adults: a systematic review of DNA methylation studies. Clin Epigenetics. 2015;7:113.
Silva CP, Kamens HM. Cigarette smoke-induced alterations in blood: a review of research on DNA methylation and gene expression. Exp Clin Psychopharmacol. 2021;29(1):116–35.
Zeilinger S, Kuhnel B, Klopp N, Baurecht H, Kleinschmidt A, Gieger C, et al. Tobacco smoking leads to extensive genome-wide changes in DNA methylation. PLoS ONE. 2013;8(5):e63812.
Tsaprouni LG, Yang TP, Bell J, Dick KJ, Kanoni S, Nisbet J, et al. Cigarette smoking reduces DNA methylation levels at multiple genomic loci but the effect is partially reversible upon cessation. Epigenetics. 2014;9(10):1382–96.
Lee KW, Pausova Z. Cigarette smoking and DNA methylation. Front Genet. 2013;4:132.
Guerrero-Preston R, Goldman LR, Brebi-Mieville P, Ili-Gangas C, Lebron C, Witter FR, et al. Global DNA hypomethylation is associated with in utero exposure to cotinine and perfluorinated alkyl compounds. Epigenetics. 2010;5(6):539–46.
Ferrari S, Karasik D. Gene-diet interactions on bone. In: Holick MF, Nieves JW, editors. Nutrition and bone health. Springer: New York; 2015. p. 21–36.
Balliu B, Carcamo-Orive I, Gloudemans MJ, Nachun DC, Durrant MG, Gazal S, et al. An integrated approach to identify environmental modulators of genetic risk factors for complex traits. Am J Hum Genet. 2021;108(10):1866–79.
Kanherkar RR, Bhatia-Dey N, Csoka AB. Epigenetics across the human lifespan. Front Cell Dev Biol. 2014;2:49.
Golbabapour S, Abdulla MA, Hajrezaei M. A concise review on epigenetic regulation: insight into molecular mechanisms. Int J Mol Sci. 2011;12(12):8661–94.
Rotondo JC, Mazziotta C, Lanzillotti C, Tognon M, Martini F. Epigenetic dysregulations in merkel cell polyomavirus-driven merkel cell carcinoma. Int J Mol Sci. 2021;22(21):11464.
Moller AMJ, Delaisse JM, Olesen JB, Canto LM, Rogatto SR, Madsen JS, et al. Fusion potential of human osteoclasts in vitro reflects age, menopause, and in vivo bone resorption levels of their donors-a possible involvement of DC-STAMP. Int J Mol Sci. 2020;21(17):1–16.
Moller AMJ, Delaisse JM, Olesen JB, Madsen JS, Canto LM, Bechmann T, et al. Aging and menopause reprogram osteoclast precursors for aggressive bone resorption. Bone Res. 2020;8:11.
Wu T, Hu E, Xu S, Chen M, Guo P, Dai Z, et al. clusterProfiler 4.0: a universal enrichment tool for interpreting omics data. Innovation (N Y). 2021;2(3):100141.
Qiu C, Yu F, Su K, Zhao Q, Zhang L, Xu C, et al. Multi-omics data integration for identifying osteoporosis biomarkers and their biological interaction and causal mechanisms. iScience. 2020;23(2):100847.
Cairoli E, Eller-Vainicher C, Ulivieri FM, Zhukouskaya VV, Palmieri S, Morelli V, et al. Factors associated with bisphosphonate treatment failure in postmenopausal women with primary osteoporosis. Osteoporos Int. 2014;25(4):1401–10.
Borgstrom F, Karlsson L, Ortsater G, Norton N, Halbout P, Cooper C, et al. Fragility fractures in Europe: burden, management and opportunities. Arch Osteoporos. 2020;15(1):59.
Kanis JA, Barton IP, Johnell O. Risedronate decreases fracture risk in patients selected solely on the basis of prior vertebral fracture. Osteoporos Int. 2005;16(5):475–82.
McCloskey E, Selby P, Davies M, Robinson J, Francis RM, Adams J, et al. Clodronate reduces vertebral fracture risk in women with postmenopausal or secondary osteoporosis: results of a double-blind, placebo-controlled 3-year study. J Bone Miner Res. 2004;19(5):728–36.
Dunford JE, Thompson K, Coxon FP, Luckman SP, Hahn FM, Poulter CD, et al. Structure-activity relationships for inhibition of farnesyl diphosphate synthase in vitro and inhibition of bone resorption in vivo by nitrogen-containing bisphosphonates. J Pharmacol Exp Ther. 2001;296(2):235–42.
Ali N, Jurczyluk J, Shay G, Tnimov Z, Alexandrov K, Munoz MA, et al. A highly sensitive prenylation assay reveals in vivo effects of bisphosphonate drug on the Rab prenylome of macrophages outside the skeleton. Small GTPases. 2015;6(4):202–11.
Liu S, Sahid MNA, Takemasa E, Maeyama K, Mogi M. Zoledronate modulates intracellular vesicle trafficking in mast cells via disturbing the interaction of myosinVa/Rab3a and sytaxin4/VAMP7. Biochem Pharmacol. 2018;151:18–25.
Sato M, Grasser W, Endo N, Akins R, Simmons H, Thompson DD, et al. Bisphosphonate action. Alendronate localization in rat bone and effects on osteoclast ultrastructure. J Clin Invest. 1991;88(6):2095–105.
Scala R, Maqoud F, Antonacci M, Dibenedetto JR, Perrone MG, Scilimati A, et al. Bisphosphonates targeting ion channels and musculoskeletal effects. Front Pharmacol. 2022;13:837534.
Yu Z, Surface LE, Park CY, Horlbeck MA, Wyant GA, Abu-Remaileh M, et al. Identification of a transporter complex responsible for the cytosolic entry of nitrogen-containing bisphosphonates. Elife. 2018;7:e36620.
Rodan GA, Reszka AA. Bisphosphonate mechanism of action. Curr Mol Med. 2002;2(6):571–7.
Bivi N, Romanello M, Harrison R, Clarke I, Hoyle DC, Moro L, et al. Identification of secondary targets of N-containing bisphosphonates in mammalian cells via parallel competition analysis of the barcoded yeast deletion collection. Genome Biol. 2009;10(9):R93.
Alakangas A, Selander K, Mulari M, Halleen J, Lehenkari P, Monkkonen J, et al. Alendronate disturbs vesicular trafficking in osteoclasts. Calcif Tissue Int. 2002;70(1):40–7.
Maunakea AK, Chepelev I, Cui K, Zhao K. Intragenic DNA methylation modulates alternative splicing by recruiting MeCP2 to promote exon recognition. Cell Res. 2013;23(11):1256–69.
Jones PA. Functions of DNA methylation: islands, start sites, gene bodies and beyond. Nat Rev Genet. 2012;13(7):484–92.
Xia Y, Liu N, Xie X, Bi G, Ba H, Li L, et al. The macrophage-specific V-ATPase subunit ATP6V0D2 restricts inflammasome activation and bacterial infection by facilitating autophagosome-lysosome fusion. Autophagy. 2019;15(6):960–75.
Luzio JP, Pryor PR, Bright NA. Lysosomes: fusion and function. Nat Rev Mol Cell Biol. 2007;8(8):622–32.
Thompson K, Rogers MJ, Coxon FP, Crockett JC. Cytosolic entry of bisphosphonate drugs requires acidification of vesicles after fluid-phase endocytosis. Mol Pharmacol. 2006;69(5):1624–32.
Sharpe LJ, Brown AJ. Controlling cholesterol synthesis beyond 3-hydroxy-3-methylglutaryl-CoA reductase (HMGCR). J Biol Chem. 2013;288(26):18707–15.
Fukami M, Horikawa R, Nagai T, Tanaka T, Naiki Y, Sato N, et al. Cytochrome P450 oxidoreductase gene mutations and Antley-Bixler syndrome with abnormal genitalia and/or impaired steroidogenesis: molecular and clinical studies in 10 patients. J Clin Endocrinol Metab. 2005;90(1):414–26.
Aryee MJ, Jaffe AE, Corrada-Bravo H, Ladd-Acosta C, Feinberg AP, Hansen KD, et al. Minfi: a flexible and comprehensive Bioconductor package for the analysis of Infinium DNA methylation microarrays. Bioinformatics. 2014;30(10):1363–9.
Fortin JP, Labbe A, Lemire M, Zanke BW, Hudson TJ, Fertig EJ, et al. Functional normalization of 450k methylation array data improves replication in large cancer studies. Genome Biol. 2014;15(12):503.
Du P, Zhang X, Huang C-C, Jafari N, Kibbe WA, Hou L, et al. Comparison of Beta-value and M-value methods for quantifying methylation levels by microarray analysis. BMC Bioinform. 2010;11(1):1–9.
Rahmani E, Zaitlen N, Baran Y, Eng C, Hu D, Galanter J, et al. Sparse PCA corrects for cell type heterogeneity in epigenome-wide association studies. Nat Methods. 2016;13(5):443–5.
Benjamini Y, Drai D, Elmer G, Kafkafi N, Golani I. Controlling the false discovery rate in behavior genetics research. Behav Brain Res. 2001;125(1–2):279–84.
Formosa MM, Bergen DJM, Gregson CL, Maurizi A, Kämpe A, Garcia-Giralt N, et al. A roadmap to gene discoveries and novel therapies in monogenic low and high bone mass disorders. Front Endocrinol. 2021;12:709711.
Rauner M, Foessl I, Formosa MM, Kague E, Prijatelj V, Lopez NA, et al. Perspective of the GEMSTONE consortium on current and future approaches to functional validation for skeletal genetic disease using cellular, molecular and animal-modeling techniques. Front Endocrinol. 2021;12:731217. https://doi.org/10.3389/fendo.2021.731217.
We wish to thank Henning Boldt for performing the DNA methylation measurements and Salah Masmoudi for setting up the database allowing us to compare the results of DNA methylation with demographic and experimental information of the corresponding participants. We acknowledge Malka Kitayner (PhD) for her analytical support. We also wish to thank Jacob Bastholm Olesen for excellent technical assistance. This work was promoted and facilitated by the membership of David Karasik and Kent Søe in the COST Action CA18139 GEMSTONE  (Genomics of MusculoSkeletal traits Translational Network), in particular their membership of Working Group 4 .
We wish to acknowledge funding from the Research Counsel of Lillebaelt Hospital; the Region of Southern Denmark (15/24819); the Department of Regional Health Research, University of Southern Denmark; the Mrs Astrid Thaysen Fund (ATL 16/02); and the Aase & Ejnar Danielsen Fund (10-001835 and 19-10-0352). The replication study was partially supported or benefited by grants from the National Institutes of Health, United States (R01AR069055, R01AG061917, U19AG055373, and P20GM109036). DK was supported by grant ISF-1121/19 from Israel Science Foundation.
Ethics approval and consent to participate
Ethical approval on the 11th of May 2015 by the Danish Regional Scientific Ethics Committee (S-20150059) and signed informed consent. The independent validation study  was approved by The Institutional Review Boards for Human Investigation at Tulane University (New Orleans, USA), and the signed informed-consent documents were obtained from all study participants.
Consent for publication
The authors declare that they have no competing interests.
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Additional file 1
. Supplementary Figure 1. Histograms of the success ratio distributions in each phenotype grouping for chosen maximal distance. X-axis, log2 of success ratio (proportion, #successful SNPs/#all SNPs): “Success” means that the distance between this SNP and a CpG in the 59 CpGs set is less than stated maximal distance. Series of incremental distances from 50,000 bp to 1000 bp are shown. Vertical axis of the histogram represents the number of SNPs in each log2 of success ratio’s bin. Data shown corresponds to only SNPs associated with: A) phenotype group 1; B) phenotype group 2; C) phenotype group 3; D) phenotype group 4.
Additional file 2
. Supplementary Figure 2. Histograms of the success ratio distributions in each phenotype grouping for chosen maximal distance. X-axis, log2 of success ratio (proportion, #successful SNPs/#all SNPs): “Success” means that the distance between this SNP and a CpG in the 59 CpGs set is less than stated maximal distance. Series of incremental distances from 50,000 bp to 1000 bp are shown. Vertical axis of the histogram represents the number of SNPs in each log2 of success ratio’s bin. Data shown corresponds to SNPs in all other groups excluding: A) phenotype group 1; B) phenotype group 2; C) phenotype group 3; D) phenotype group 4.
Additional file 3
. Supplementary Table 1. Output of EWAS statistical results (p<0.05) arranged in the order from left to the right for the interaction effect, main effect of IC50 and smoking, followed by annotations for each CpGs in the table.
Additional file 4
. Supplementary Table 2. Validation of top 59 significant smoking-dependent CpGs in an independent multi-omics BMD study . Significant q-values are highlighted in bold, while highly significant q-values are also highlighted in red.
Additional file 5
. Supplementary Table 3. Validation of top CpGs showing main effect on IC50 in an multi-omics independent BMD study . Significant q-values are highlighted in bold, while highly significant q-values are also highlighted in red.
Additional file 6
. Supplementary Table 4. The number and list of phenotypes in each phenotype group. Each sheet is labelled according to the corresponding “phen#” group.
Additional file 7
. Supplementary Table 5. This table shows the raw data for each CpG and each Phen#, success ratio was calculated as a proportion of #successful SNPs for that CpG to #all SNPs associated with that phenotype group, Phen#. “Success” means that the distance between location of a SNP and location of the CpG is less than the stated max distance. Series of incremental distances from 50,000 bp to 1,000 bp were tested, and then compared across the phenotype groups. The one-tailed Mann-Whitney U test was used to compare differences in success ratios between Phen1 and Phen 2 and “negative control” Phen4 as shown in Table 4.
Rights and permissions
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated in a credit line to the data.
About this article
Cite this article
Tan, Q., Møller, A.M.J., Qiu, C. et al. A variability in response of osteoclasts to zoledronic acid is mediated by smoking-associated modification in the DNA methylome. Clin Epigenet 15, 42 (2023). https://doi.org/10.1186/s13148-023-01449-1
- Zoledronic acid
- DNA methylation
- Association studies