Pharmacogenomic Effects of β-Blocker Use on Femoral Neck Bone Mineral Density.

CONTEXT: Recent studies have shown that β-blocker (BB) users have a decreased risk of fracture and higher bone mineral density (BMD) compared to nonusers, likely due to the suppression of adrenergic signaling in osteoblasts, leading to increased BMD. There is also variability in the effect size of BB use on BMD in humans, which may be due to pharmacogenomic effects.
OBJECTIVE: To investigate potential single-nucleotide variations (SNVs) associated with the effect of BB use on femoral neck BMD, we performed a cross-sectional analysis using clinical data, dual-energy x-ray absorptiometry, and genetic data from the Framingham Heart Study's (FHS) Offspring Cohort. We then sought to validate our top 4 genetic findings using data from the Rotterdam Study, the BPROOF Study, the Malta Osteoporosis Fracture Study (MOFS), and the Hertfordshire Cohort Study.
METHODS: We used sex-stratified linear mixed models to determine SNVs that had a significant interaction effect with BB use on femoral neck (FN) BMD across 11 gene regions. We also evaluated the association of our top SNVs from the FHS with microRNA (miRNA) expression in blood and identified potential miRNA-mediated mechanisms by which these SNVs may affect FN BMD.
RESULTS: One variation (rs11124190 in HDAC4) was validated in females using data from the Rotterdam Study, while another (rs12414657 in ADRB1) was validated in females using data from the MOFS. We performed an exploratory meta-analysis of all 5 studies for these variations, which further validated our findings.
CONCLUSION: This analysis provides a starting point for investigating the pharmacogenomic effects of BB use on BMD measures.

Osteoporosis is a skeletal condition that causes bones to become fragile, resulting in an increased risk of fracture and decreased bone mineral density (BMD). This disorder affects more than 10 million individuals in the United States and results in more than 2 million osteoporotic fractures per year, with the annual hospital cost of osteoporotic fractures exceeding $28 billion dollars [1-4]. Several studies have found an association between β-blocker (BB) use, decreased risk of fracture, and higher BMD [5-8], including a pilot randomized trial [9] This association is thought to be mediated, at least in part, by attenuation of adrenergic signaling in osteoblasts (Fig. 1) [10-13]. In particular, it has been found that norepinephrine signaling activates β-adrenergic receptors in osteoblasts leading to signaling through cyclic adenosine 5′-monophosphate and protein kinase A, resulting in the activation of ATF4. ATF4 is a transcription factor that induces transcription of TNFSF11 (RANKL). TNFSF11 (RANKL) is secreted from osteoblasts and binds to either TNFRSF11A (RANK), a receptor on the surface of osteoclasts, or osteoprotegerin (OPG), a soluble decoy receptor produced by osteoblasts. TNFSF11 (RANKL) signaling through TNFRSF11A (RANK) leads to increased bone resorption due to increased osteoclast activity and differentiation.

Figure 1.

Adrenergic signaling in bone. Norepinephrine (NE) binds to β-adrenergic receptors, stimulating adrenergic signaling through 3′,5′-cyclic adenosine 5′-monophosphate (cAMP) and PKA. This results in the activation of ATF4, a transcription factor that triggers the transcription of TNFSF11 (RANKL). HDAC4 is a histone deacetylase that further acts to stabilize ATF4. TNFSF11 (RANKL) is secreted by osteoblasts and binds to TNFRSF11A (RANK) receptors on osteoclasts or osteoprotegerin (OPG)-soluble decoy receptors. Activation of TNFRSF11A (RANK) then stimulates osteoclast differentiation, leading to bone resorption. β-Blockers competitively bind to β-adrenergic receptors, blocking signaling by norepinephrine.

While many studies have found an association between BB use and bone outcomes, there is variability in the effect size [5-7], and some negative studies [14-16]. We hypothesize that genetic variation may contribute to this variability given the large genetic component of BMD itself [17], and given the pharmacogenetic effects found for cardiovascular outcomes [18-21], with demonstrated associations in the β-adrenergic receptor genes [12, 18, 22-25]. However, associations between these single-nucleotide variation (SNVs; formerly single-nucleotide polymorphisms [SNPs]) and BMD have not been demonstrated, as a recent analysis by Veldhuis-Vlug et al showed that nonsynonymous SNVs in ADRB2 were not significantly associated with BMD or fracture risk [26]. Previous genome-wide association studies (GWAS) of BMD and osteoporosis have been performed [27-29], but there have not been previous pharmacogenomic studies evaluating the effect of BB use on BMD, although other studies have found SNVs associated with BMD and osteoporosis to map to or near genes involved in adrenergic signaling [10, 27, 30].

In addition to genetic variations, we have previously found circulating microRNAs (miRNAs) to be associated with BB use and BMD and revealed potential miRNA-mediated mechanisms by which BB use influences BMD, including attenuation of adrenergic signaling in osteoblasts [31]. miRNAs are small (~ 22 nucleotides), noncoding RNAs that act on target messenger RNAs (mRNAs) to inhibit protein expression through mRNA degradation and translational inhibition [32]. Circulating miRNAs have been used to develop hypotheses regarding underlying mechanisms in many applications, including cardiovascular disease (CVD) and cancer etiology, variation in handgrip strength, and response to antidepressant treatment [33-35]. Several circulating miRNAs have been implicated as potential biomarkers of osteoporosis and BB treatment response [36-40], and miRNAs have also been assessed in association with GWAS signals using expression quantitative loci analysis in many outcomes to discover potential mechanisms and biomarkers of these conditions [41-43].

To discover genetic variations associated with the effect of BB use on femoral neck (FN) BMD, we sought to evaluate genetic variations that map to or near genes involved in adrenergic signaling in bone and that interact with BB use in their association with FN BMD. These candidate genes were chosen based on what is currently known about the effect of adrenergic signaling on osteoblasts in bone as has been previously described (see Fig. 1) [10-13]. FN BMD was chosen as an outcome variable because of its clinical importance and its use to evaluate fracture risk of patients using the World Health Organization fracture risk assessment (FRAX) tool [44]. We also sought to determine putative underlying miRNA-mediated pathways involved in this association.

We hypothesize that genetic variations in these candidate genes can partly explain the variation in FN BMD among BB users. Furthermore, we hypothesize that genetic variations may affect miRNA-mediated mechanisms underlying the association between BB use and FN BMD. To test these hypotheses, we followed the analytic plan shown in Supplementary Figure S1 (Supplementary Material, http://clary.mmcri.cloud/JES_Nevola_et_all_supplement.pdf). We used linear mixed modeling followed by conditional joint analysis (COJO) to analyze this genetic association in these candidate genes using clinical data from the Framingham Heart Study (FHS). Since β-adrenergic signaling and the effect of BB use has previously been shown to have sex-specific effects [45-49], and also BMD and osteoporosis prevalence vary by sex [50-53], we used a sex-stratified model for our genetic association studies. We submitted 4 SNVs for validation and obtained validation for 2 in independent studies, and additionally performed exploratory meta-analyses across cohorts. To generate functional hypotheses, we analyzed individual miRNAs that were associated with these variations and BB use to identify candidate mechanisms that were altered in the presence of the alternative alleles.

Materials and Methods

Study Sample

Data for this cross-sectional analysis were made available from dbGaP through approved request number 1302685-1 [54]. The FHS is an ongoing, 3-generation, community-based study. For this study, we focused on members of the Offspring Cohort, which includes the children of the original cohort and their spouses. At each FHS examination, age, height, body mass index (BMI), and extensive questionnaires were obtained according to standardized protocols. Most of the members of the Offspring Cohort were enrolled in the ancillary Framingham Osteoporosis Study in 2002 [55]. BMD was measured at the hip (FN, trochanter, and total femur) and lumbar spine (average BMD of L2-L4) in grams divided by centimeters squared (g/cm2) using a GE Lunar Prodigy dual-energy x-ray absorptiometer. For this analysis, 1527 individuals were included based on being a member of the Framingham Offspring Cohort who attended examination cycle 8 (2005-2008, n = 3021), having BMD data that were assessed after the exam 8 date when BB use was assessed, and having genetic data available. Genetic data were collected and imputed as previously described [56, 57]. In brief, genotypes were measured using the Affymetrix 500K and 50K Human Gene Focused Panels. Genetic variations’ positions were based on the GRCh37/hg19 assembly from February 2009. Imputation was based on the Haplotype Reference Consortium reference panel release 1. The panel included only autosomes with 39 235 157 sites, of which 39 210 718 sites were included in the data set returned by the Michigan Imputation Server with high-quality imputation. Multiallelic sites were excluded from our analysis. The imputed SNVs’ value ranged from 0 to 2, referring to the predicted dosage of the alternative allele.

Medication Assessment

Medication usage, including oral BB use as the primary exposure and other medications related to bone in sensitivity analyses, was measured using a medication questionnaire in which the medication name, strength, route, and frequency (day/week/month/year) were recorded by directly viewing the medication bottle during the exam 8 (2005-2008) visit, excluding as-needed use. We categorized BB users as β1-selective for the chemical group “Beta blocking agents, selective” and as “non-selective” for the chemical groups “Beta blocking agents, non-selective” or “Alpha and beta blocking agents.” We computed BB daily dose for each patient and for each drug by converting the strength and frequency to a daily dose. We divided this calculated daily dose by the World Health Organization–determined defined daily dose [58] to get a standardized dose in units of defined daily dose for that drug.

We additionally recorded use of other bone-related medications for sensitivity analyses examining confounding by these variables. Medications for therapy group “bone diseases” or with the chemical name “raloxifene” were considered bone disease drugs. The chemical names for the therapy group bone diseases consisted mostly of alendronic acid and risedronic acid, with etidronic acid, ibandronic acid, ipriflavone, pamidronic acid, and zoledronic acid also included. We also noted oral steroid use (including chemical groups “corticosteroids” and “glucocorticoids”). Treatment for hypertension, lipids, or diabetes was recorded as part of the exam 8 visit. Prior CVD was determined from an adjudicated file of cardiovascular events recorded prior to exam 8.

All participants provided informed consent, and the examination protocols were approved by the Boston University Medical Center Institutional Review Board and the Hebrew SeniorLife Institutional Review Board.

MicroRNA Expression Profiling

Whole blood from fasting morning samples was used for miRNA profiling, which was obtained at exam 8, coincident with the BB use ascertainment, and just before BMD measurement and stored at –80 °C. Several studies have used these miRNA data in association with BMD, BB use, and other phenotypes [31, 35, 59, 60]. In brief, the high-throughput Gene Expression Core Laboratory at the University of Massachusetts Medical School profiled commercially available TaqMan miRNA assays [35, 59, 60]. A subset of the 754 miRNAs profiled in 600 FHS participants was further profiled in additional FHS Offspring Cohort members using quantitative real-time polymerase chain reaction (RT-qPCR). A total of 333 miRNAs had a measurable cycle threshold (Ct) value in at least 5% of participants. A higher Ct value reflects a lower miRNA expression value. The FHS Systems Approach to Biomarker Research in Cardiovascular Disease Initiative Steering Committee previously reviewed all quality control measures and noted that more than 95% of the data points had coefficients of variation of less than 10% (mean ~ 4%) [60]. Of the 1527 individuals who were included in the genetic association study, 1304 had miRNA data available for miRNA association analysis. We modeled technical sources of variation in miRNA concentration (Ct) values including RNA quality, RNA concentration, and 260/280 ratio (ratio of absorbance at 260 and 280 nm using a spectrophotometer) as previously described [35]. Briefly, we categorized each technical variable by decile and included them as factor variables in our models to account for nonlinear effects.

Identifying Genes of Interest

For our candidate gene list, we selected genes involved in adrenergic signaling in osteoblasts as detailed by Elefteriou et al [10]. We added TNFRSF11A (RANK) and OPG to this list of genes because TNFSF11 (RANKL) binds to the receptors encoded by these genes. Our pathway, therefore, starts at β-adrenergic receptors and ends at TNFSF11 (RANKL) receptors TNFRSF11A (RANK) and OPG (see Fig. 1). To further filter our list, we also required that at least one SNV that mapped to or near each gene have a suggestive association with eBMD (BMD estimated by quantitative ultrasound of the heel), FN BMD, or lumbar spine BMD in a previous GWAS studies as reported in the Musculoskeletal Genomics Knowledge Portal [30]. A suggestive association was defined as having at least 1 SNV within the coding region ± 100 kb that is associated with the phenotype with a P value less than 5 × 10–4, parameters that have been set by the Musculoskeletal Genomics Knowledge Portal. Similar parameters have been used as suggestive P values in previous studies, ranging from 1 × 10–4 to 1 × 10–6 [61, 62]. Our final gene list contained 13 genes: 3 β-adrenergic receptor genes (ADRB1, ADRB2, ADRB3), 5 PKA subunits (PRKACB, PRKAR1A, PRKAR1B, PRKAR2A, PRKAR2B), HDAC4, ATF4, TNFSF11 (RANKL), TNFRSF11A (RANK), and OPG.

Identifying Single-Nucleotide Variations of Interest

We performed our analysis in a 2-stage design, in which we used Framingham as the discovery cohort and then sought to validate our top SNVs in replications cohorts. This strategy has been used in previous pharmacogenetic studies [63]. Variations had previously been excluded if they satisfied any of the following criteria: Hardy-Weinberg equilibrium value P value less than 1 × 10– 6, call rate less than 96.9%, minor allele frequency less than 0.01, number of mendelian errors greater or equal to 1000, or at locations that did not map to GRCh37 [56, 57]. Well-imputed SNVs were determined across the genome by filtering for an R2 value greater than 0.8 as provided by the Haplotype Reference Consortium after imputation. Then SNVs from 2 kB upstream of the gene region to 0.5 kB downstream of the gene region were extracted for further analysis according to RefSeq [64-67], filtering out polyallelic SNVs. These SNVs were then filtered for having a minor allele frequency greater than 0.05 in our population of 1527 individuals. This resulted in 1482 SNVs across 11 genes. ADRB3 and PRKAR1B did not have any SNVs that met our filtering criteria.

Modeling Interaction Effect Between Genotype and β-Blocker Use on Femoral Neck Bone Mineral Density

The 1482 variations were analyzed in 1527 individuals for an association with FN BMD that was modified by BB use using a linear mixed model, stratifying based on sex, and adjusting for interrelatedness between individuals by modeling a kinship matrix as a random effect (lmekin function in coxme package in R). We performed a sex-stratified analysis since β-adrenergic signaling and the effect of BB use has previously been shown to have sex-specific effects [45-49], and also BMD and osteoporosis prevalence vary by sex [50-53]. The female-only model adjusted for age, height, BMI, and current estrogen use, and the male-only regression model adjusted for age, height, and BMI. Menopausal status was not adjusted for because more than 99.6% of our female cohort was postmenopausal; only 3 women were premenopausal (0.37%). FN BMD measurement was used as the dependent variable, while allele dosage, BB use, and their interaction for each SNV were modeled, and the interaction effect estimate and P value were used as the parameters of interest. FN BMD was chosen as the dependent variable because of its clinical importance, especially in the calculation of the FRAX score, and FN BMD was available in more participants than other BMD sites, such as lumbar spine BMD [44]. We also focused on a single skeletal site because of the limited power for detecting interaction in a pharmacogenetic study, and the multiple testing penalty that would ensue with multiple phenotypes.

To account for linkage disequilibrium (LD) between SNVs, we performed a COJO using GCTA (GCTA-COJO) [68, 69], and filtered for a P value of less than .05 after the COJO. GCTA-COJO was used to perform a stepwise model selection procedure to select independently associated SNVs. Default parameters were used with the exception of threshold P value, which was set to .05, and difference in allele frequency between summary statistics and LD reference sample, which was set to 1. Genetic data from the FHS (n = 1527) were used both as the reference sample to estimate LD as well as the data set to create the summary statistics file. We did not perform additional multiple testing adjustment in the discovery phase, and instead performed a Bonferroni correction for the number of SNVs in the validation phase, as in Singh et al [63]. Since we adjusted for covariates in our models that may have genetic components (height, BMI), we examined previous SNV associations with these covariates in prior GWAS studies to identify potential collider bias.

Determining Single-Nucleotide Variations in High Linkage Disequilibrium

We used HaploReg [70], Search Candidate cis-Regulatory Elements by ENCODE (SCREEN) [71], and LDlink [72] to explore LD among SNVs and annotations including chromatin state, previous expression quantitative loci signals, and proteins bound in ChIP-Seq experiments for our top 11 SNVs and SNVs in high LD with those SNVs (r2 = 0.8) as calculated using the European population of the 1000 Genomes Project using HaploReg or LDLink. Correlation between SNVs was also calculated within our cohort using the R function cor.

MicroRNA Association Analysis

The 11 genetic variations that were found to be significant in our pharmacogenomic association model were analyzed in 1304 of the 1527 individuals who had available miRNA data. We determined association between SNVs and 333 miRNAs assayed using qPCR data. Associated miRNAs were determined using a linear mixed model, stratifying based on sex and adjusting for interrelatedness between individuals by modeling a kinship matrix as a random effect (lmekin function in coxme package in R [73]). The female-only model adjusted for age, height, BMI, current estrogen use, and miRNA technical variables (RNA concentration, RNA quality, and RNA 260/280 ratio, a measure of purity of the RNA), and the male-only regression model adjusted for age, height, BMI, and the same miRNA technical variables. Isolation batch effect was not included as a covariate because of power restraints. miRNA expression level as measured by qPCR was used as the dependent variable, while allele dosage, BB use, and their interaction were modeled for each SNV.

For many participants, some miRNAs were not expressed at a detectable level. Therefore, for each miRNA in each participant, the expression level was redefined as a discrete variable, X = 1 if undetectable and X = 0 if detectable. For miRNA expressed in more than 5% but less than 10% of samples, we modeled the discrete expression value, and for miRNA expressed in more than 90% of samples, we modeled the continuous expression value. For miRNA expressed in 10% to 90% of samples, the discrete model and continuous expression model were both computed with the final P value determined using the Fisher method.

MicroRNA Target Determination

mRNA targets of significant miRNA were determined using the get_multimiR function in the multimiR R Package [74], which queried miRNA-target databases to determine validated targets of each miRNA (miRecord, miRTarBase, TarBase).

Validation Analysis

Our top 4 SNVs, 2 in each sex, were submitted for validation in 4 independent population-based cohorts: the Rotterdam Study, the BPROOF study, the Malta Osteoporosis Fracture Study (MOFS), and the Hertfordshire Cohort Study. The variations were assessed for a significant interaction effect with BB use, adjusting for multiple testing using prespecified thresholds, using linear regression modeling in validation cohorts. The female-only model adjusted for age, height, BMI, and current estrogen use, and the male-only regression model adjusted for age, height, and BMI. A summary of validation cohorts’ data and methodology can be found in Supplementary Table S1 (Supplementary Material, http://clary.mmcri.cloud/JES_Nevola_et_all_supplement.pdf). Analysis of the Rotterdam study also adjusted for cohort effect in all models. The prespecified criteria for an SNV to be validated was P less than .05/4 in cohorts with both sexes and P less than .05/2 in cohorts with only 1 sex.

Meta-Analysis Methods

Meta-analysis was performed across the 5 cohorts using the metagen function in the meta package [75]. A fixed-effect and random-effect meta-analysis was performed based on the effect estimates and their SEs. The inverse variance method was used for pooling. Forest plots were generated using the forest function in the meta package [75]. The fixed-effect model estimates were used unless significant heterogeneity, as calculated by the I2 statistic, was observed (P value of I2 statistic < 0.05).

Results

Eleven Single-Nucleotide Variations Found to Be Significant in 6 Genes in Discovery Sample

Characteristics of the study cohort are given in Supplementary Table S2 (http://clary.mmcri.cloud/JES_Nevola_et_all_supplement.pdf), including estrogen usage rate, which was 8.9% in women, and use of medication for bone disease, which was 21.3% for women and 2.8% for men. We analyzed 1482 SNVs across 13 genes related to adrenergic signaling in bone using genetic data from the FHS (Supplementary Table S2, http://clary.mmcri.cloud/JES_Nevola_et_all_supplement.pdf). Eleven SNVs in 6 genes were found to have a significant interaction effect with BB use on FN BMD (P < .05) after GCTA-COJO analysis (Table 1). Five SNVs were found to be significant in the female-only model, while the other 6 were significant in the male-only model. There was no overlap in significant SNVs across sexes. Most of these SNVs were intronic variations, except for rs12414657 (ADRB1), which is an upstream transcript variation, and rs13393217 (TNFSF11 or RANKL), which is a 3 prime untranslated region variation. We also looked at the functional annotation of highly correlated SNVs using HaploReg, SCREEN, and LDlink and by performing correlation analysis within the Framingham cohort. Of these 11 SNVs, only one had a nonsynonymous SNV in high LD, rs12414657 (ADRB1), which is highly correlated with rs1801252, a missense SNV that codes for a serine to glycine shift at the 49th amino acid in ADRB1.

Table 1.

Significant single-nucleotide variations using conditional joint analysis using GCTA

Gene (Ref/Alt)	Position hg19	rsID	Model	Effect size from conditional analysis	SE from conditional analysis	P from conditional analysis
ADRB1 (T/C)	115803375	rs12414657	Female	0.0431	0.0216	.05
HDAC4 (C/A)	239972561	rs13393217	Female	0.0440	0.0194	.02
HDAC4 (C/G)	240223080	rs11124190	Female	0.0489	0.0188	.009
HDAC4 (G/A)	240050108	rs145900122	Male	0.0872	0.0343	.01
HDAC4 (G/A)	240112014	rs3791554	Male	–0.0576	0.0257	.02
PRKACB (A/G)	84682179	rs970318	Male	0.0373	0.0162	.02
PRKAR2B (G/A)	106736732	rs6952920	Female	0.0401	0.0136	.003
TNFRSF11A (RANK) (T/C)	60025809	rs72933609	Female	0.0703	0.0293	.02
TNFRSF11A (RANK) (G/A)	60001153	rs34170507	Male	0.0695	0.0170	4.4 × 10 –5
TNFRSF11A (RANK) (C/T)	60026732	rs6567268	Male	0.0484	0.0173	.005
TNFSF11 (RANKL) (T/C)	43177169	rs9533166	Male	–0.0310	0.0158	.05

SNVs that met a P value less than .05 cutoff using GCTA-COJO analysis, including the gene the SNV is located in or near and the reference and alternative alleles, the position of the SNV in the hg19 Genome Build, the rsID of the SNV, and the model in which the SNV was significant (female-only or male-only model). The effect size, SE, and P value were determined using conditional joint analysis using the summary statistics from the linear mixed model analysis. The linear mixed model included the interaction effect between the alternative allele dosage of the SNV and β-blocker use and its effect on femoral neck bone mineral density, adjusting for covariates and modeling interrelatedness between individuals using a kinship matrix. The summary statistics for all SNVs were then used to perform GCTA-COJO analysis. SNVs chosen for validation are included in bold.

Abbreviations: GCTA-COJO, conditional joint analysis using GCTA; SNV, single-nucleotide variation (formerly single-nucleotide polymorphism [SNP]).

Single-Nucleotide Variations for Validation

We chose to validate 4 of the 11 SNVs in an external cohort, 2 in females and 2 in males, to limit our multiple testing burden, which was strictly controlled in our validation cohorts. These SNVs are indicated in bold in Table 1 and Supplementary Table S3. Our rationale for validation is detailed in Supplementary Figure S2 (http://clary.mmcri.cloud/JES_Nevola_et_all_supplement.pdf). Of the 11 SNVs with a P value less than .05 after conditional analysis, none were nonsynonymous SNVs, but rs12414657 (ADRB1) was in high LD with a nonsynonymous SNV, so this SNV was chosen for validation in females. The SNVs that mapped to or near PKA subunit genes (rs970318 and rs6952920) were excluded from validation because PKA is involved in many different processes, and the SNVs could not be mapped to a role in β-adrenergic signaling or BMD. The most significant SNVs in each sex were then chosen to reach 2 SNVs per sex. These SNVs were rs11124190 (HDAC4) in females and rs34170507 and rs6567268 (both in TNFRSF11A or RANK) in males. SNVs were considered validated in an external cohort if they met the following prespecified significance thresholds: P less than .0125 (0.05/4) in cohorts with both sexes and P less than .025 (0.05/2) in cohorts with only one sex to account for multiple testing. The effect estimate for the interaction effect of the alternative allele was positive for all 4 SNVs in the discovery sample (FHS), indicating higher BMD in BB users with more copies of the alternative allele compared with non-BB users. Of these 4 SNV, rs11124190 (HDAC4) had a significant interaction effect with BB use on BMD in females from the Rotterdam Study (estimate = 0.024, SE = 0.009, P = .010) (Fig. 2, Supplementary Table S3; http://clary.mmcri.cloud/JES_Nevola_et_all_supplement.pdf), and rs12414657 (ADRB1) had a significant interaction effect with BB use on BMD in females from the MOFS (estimate = 0.0576, SE = 0.0219, P = .0085) (Fig. 3, Supplementary Table S3; http://clary.mmcri.cloud/JES_Nevola_et_all_supplement.pdf). The other 2 SNVs were not significant in males in the Rotterdam, BPROOF, or Hertfordshire Cohort studies (see Supplementary Table S3; http://clary.mmcri.cloud/JES_Nevola_et_all_supplement.pdf). We also performed a meta-analysis for rs11124190 (HDAC4) and for rs12414657 (ADRB1) in all 5 studies, with a significant interaction in the fixed-effect model for both SNVs (fixed-effect model for rs11124190: estimate = 0.0166, CI, 0.0035-0.0296, P = .0128; fixed-effect model for rs12414657: estimate = 0.0168, CI, 0.0015-0.0320, P = .0314) (Figs. 2 and 3). There was no evidence of significant heterogeneity at either locus across these studies.

Figure 2.

Forest plot of meta-analysis for rs11124190 (HDAC4) in females. Meta-analysis between the Framingham Heart Study (FHS), the Rotterdam Study, the BPROOF Study, the Malta Osteoporosis Fracture Study (MOFS), and the Hertfordshire Cohort Study for rs11124190 (HDAC4) in female-only models. TE is the treatment estimate and refers to the estimate of each model; seTE refers to the SE of the treatment estimate. The weight (fixed) and weight (random) columns refer to the weighting for the fixed-effects model and the random-effects model, respectively.

Figure 3.

Forest plot of meta-analysis for rs12414657 (ADRB1) in females. Meta-analysis between the Framingham Heart Study (FHS), the Rotterdam Study, the BPROOF Study, the Malta Osteoporosis Fracture Study (MOFS), and the Hertfordshire Cohort Study for rs12414657 (ADRB1) in female-only models. TE is the treatment estimate and refers to the estimate of each model; seTE refers to the SE of the treatment estimate. The weight (fixed) and weight (random) columns refer to the weighting for the fixed-effects model and the random-effects model, respectively.

MicroRNAs Associated With Top Single-Nucleotide Variations

To determine potential miRNA-related mechanisms for these SNVs, we determined significantly associated miRNA (P < .05) with each of the top 11 SNVs. We then determined if the associated miRNAs had been previously associated with osteoporosis or BMD measures, which we term “bone-related miRNAs” (Supplementary Table S4; http://clary.mmcri.cloud/JES_Nevola_et_all_supplement.pdf) [36-38, 76-93]. We also noted the association of SNVs with miR-19a-3p and miR-186-5p because we have previously found these to be associated with BB use and BMD [31]. Finally, we determined if any of the significant miRNAs targeted the gene in which the associated SNV is located in or nearby. Of note, 8 of the 11 SNVs were associated with bone-related miRNAs, and 5 SNVs were associated with miRNAs that targeted the gene where the SNV is located (Supplementary Table S5; http://clary.mmcri.cloud/JES_Nevola_et_all_supplement.pdf).

Bone-Related MicroRNAs Associated With Top Single-Nucleotide Variations

To develop hypotheses regarding the mechanism by which these SNVs interact with BB use to influence BMD, we evaluated miRNA associated with our top 4 SNVs that we tested for validation. These SNVs were associated with at least one bone-related miRNA or a miRNA that targeted the gene in which the SNV is located (Supplementary Table S5; http://clary.mmcri.cloud/JES_Nevola_et_all_supplement.pdf). Rs12414657 (ADRB1) was associated with increased miR-19a-3p expression in female BB users (Supplementary Fig. S3A and S3B; http://clary.mmcri.cloud/JES_Nevola_et_all_supplement.pdf). We have previously found miR-19a-3p to be positively associated with BB use, total femur BMD, and lumbar spine BMD [31]. ADRB1 is also a validated target of miR-19a-3p [94]. The rs11124190 (HDAC4) variation was associated with decreased expression of miR-17-5p in female BB users (Supplementary Fig. S3C and S3D; http://clary.mmcri.cloud/JES_Nevola_et_all_supplement.pdf). miR-17-5p is associated with osteoporosis (Supplementary Table S4; http://clary.mmcri.cloud/JES_Nevola_et_all_supplement.pdf) and is a biomarker of osteoporosis and suppresses osteogenic differentiation [76]. In TNFRSF11A (RANK), rs34170507 was associated with decreased expression of miR-31-5p in male BB users (Supplementary Fig. S4A and S4B; http://clary.mmcri.cloud/JES_Nevola_et_all_supplement.pdf) and miR-31-5p suppresses osteogenic differentiation [95]. Finally, rs6567268 (TNFRSF11A or RANK) was associated with increased expression of let-7g-5p and miR-374a-5p in male BB users (Supplementary Fig. S4C, S4D, and S4E; http://clary.mmcri.cloud/JES_Nevola_et_all_supplement.pdf). Let-7g-5p and miR-374a-5p target TNFRSF11A (RANK) mRNA and suppress its expression (found using multiMir R package [74]). These miRNAs may provide insights into potential mechanisms by which BB users with the alternative allele of these genetic polymorphisms tend to have higher BMD (Fig. 4).

Figure 4.

Hypothesized microRNA (miRNA)-mediated mechanisms underlying the association between top single-nucleotide variations and bone mineral density (BMD) in β-blocker (BB) users. Female BB users with the alternative allele of rs12414657 (ADRB1) have higher expression of miR-19a-3p and higher BMD. miR-19a-3p inhibits gene targets involved in adrenergic signaling, including ADRB1 and HDAC4. This inhibition of adrenergic signaling in bone would then lead to increased BMD. Female BB users with the alternative allele of rs11124190 (HDAC4) have lower expression of miR-17-5p and higher BMD. miR-17-5p inhibits osteogenic differentiation, therefore lower expression of miR-17-5p would lead to higher BMD. Male BB users with the alternative allele for rs34170507 (TNFRSF11A [RANK]) have lower expression of miR-31-5p and higher BMD. miR-31-5p inhibits osteogenic differentiation, so lower expression of miR-31-5p should lead to higher BMD. Male BB users with the alternative allele for rs6567268 (TNFRSF11A or RANK) have higher expression of let-7g-5p and miR-374a-5p and higher BMD. Let-7g-5p and miR-374a-5p both inhibit TNFRSF11A (RANK) expression. The lower TNFRSF11A (RANK) expression would decrease bone resorption, leading to higher BMD.

Sensitivity Analyses

Since BBs are used for several treatment indications, and BB users may be taking other medications or have comorbidities that may influence BMD, we performed a series of sensitivity analyses to address potential confounding in our top 2 validated genetic variants in females. The number of individuals in each medication or comorbidity category is summarized in Supplementary Table S2 (http://clary.mmcri.cloud/JES_Nevola_et_all_supplement.pdf). First, we repeated our primary analysis excluding medications taken for bone disease (see “Materials and Methods”) and found the interaction effect of BB use and SNV (number of alternative alleles) to be almost identical to the original model (Supplementary Table S6; http://clary.mmcri.cloud/JES_Nevola_et_all_supplement.pdf) with an effect size of 0.042 (SE 0.024) for the ADRB1 SNV and 0.053 (SE 0.022) for the HDAC4 SNV, although the P value did become nonsignificant because of the loss of power for ADRB1 (P = .08) but remained significant for HDAC4 (P = .017). We additionally excluded oral steroid use and found a similar result. Next, to account for hypertension or lipid treatment, we chose to analyze our interaction models within each treatment category after excluding those treated for bone disease. The BB by SNV interaction effect estimates are shown in Supplementary Table S6; http://clary.mmcri.cloud/JES_Nevola_et_all_supplement.pdf) and show remarkable consistency with the original estimates. They are all positive, and while the P values do increase because of the reduction in subset sample sizes, they remain significant or suggestive in most cases. We also adjusted for treatment for diabetes and found a nearly identical effect estimate and found similar effect estimates when excluding those with prior CVD (see Supplementary Table S6; http://clary.mmcri.cloud/JES_Nevola_et_all_supplement.pdf).

Next, because β1-selectivity and dosage of the BB used may influence the genetic interaction, we fit models in which we compared β1-selective BB use vs no BB use and then standardized daily dose (see “Materials and Methods”) as a linear term or categorized at the median into “low” or “high” values. The β1-selective users showed slightly reduced effect sizes, and the dose model showed large and highly significant effects at both loci. Furthermore, when stratifying into low- and high-dose BB users, the interaction effect was found to be much larger and more significant in the high-dose groups for both loci. We also looked at the effects of BB use and the top 2 SNVs on total hip BMD and lumbar spine BMD (see Supplementary Table S6; http://clary.mmcri.cloud/JES_Nevola_et_all_supplement.pdf). We found the effect sizes to be similar in total hip compared with femoral neck though slightly reduced in size, although still significant in the case of the HDAC4 locus. In the case of lumbar spine, the effect sizes were slightly increased at both loci although not quite significant.

Finally, as hemolysis may be a confounder of miRNA differential expression in blood, we also performed sensitivity analyses in which we additionally adjusted for miR-451a expression, an miRNA that is associated with hemolysis, for the miRNA relevant to our top SNV candidates (miR-19a-3p for the ADRB1 SNV and miR-17-5p for the HDAC4 SNV). We found that this adjustment causes a slight decrease in effect size for each of these models (see Supplementary Table S6; http://clary.mmcri.cloud/JES_Nevola_et_all_supplement.pdf), but that the P values remain quite low and are nearly significant. We additionally tested the association between miR-451a expression and BB use, FN BMD, and miR-19a-3p and miR-17-5p, and did not find any of the associations to be significant (P > .05). Therefore, we feel it is unlikely that hemolysis has confounded our miRNA genetic association results.

Discussion

We have identified 4 SNVs in the FHS discovery cohort that show a significantly higher BMD for BB users with more copies of the alternative allele as compared to nonusers, 2 of which were validated in external cohorts (rs11124190 [HDAC4] and rs12414657 [ADRB1]) in females, and 2 of which were not validated (rs34170507 and rs6567268 in TNFRSF11A or RANK in males). This focused our genetic analysis on 3 genes involved in the adrenergic signaling pathway in bone: ADRB1, HDAC4, and TNFRSF11A (RANK). We have also identified 5 miRNAs that are associated with the interaction effect between these SNVs and BB use that are associated with genes in the adrenergic signaling pathway in bone or have previously been shown to be associated with osteoporosis. As such, we have presented 4 putative mechanisms by which these SNVs interact with BB use to influence BMD (see Fig. 4).

Of these 4 SNVs, rs11124190 (HDAC4) was validated in the Rotterdam Study and rs12414657 (ADRB1) was validated in the MOFS. The minor allele frequency of rs11124190 (HDAC4) in European cohorts is 0.15, while the minor allele frequency of rs12414657 (ADRB1) is 0.14 [70]. The effect size of the HDAC4 SNV is 0.048 g/cm2 and for the ADRB1 SNV is 0.043 g/cm2, which represent the difference in FN BMD between BB users and nonusers for each additional copy of the alternative allele. A magnitude of 0.043 g/cm2 represents a 4.7% difference in BMD for females who have an average of 0.91 g/cm2 (see Supplementary Table S2; http://clary.mmcri.cloud/JES_Nevola_et_all_supplement.pdf) in this study. Considering that the average annual loss of FN BMD is 0.6% in older men and women (average age 75 years) [96], this effect size represents close to 8 years of BMD loss due to aging, although rates vary by age, sex, and other factors, thus these effect sizes are clinically significant.

The intronic HDAC4 SNV, rs11124190, is a variation in high LD with other intronic variations in HDAC4. According to HaploReg and SCREEN, this SNV has not been reported to be associated with methylation or acetylation histone modification in osteoblast primary cells, nor did it overlap with any cis-regulatory elements in other cells [70, 71]. This SNV has not been previously reported in association with other traits. Rs12414657 is a variation 430 bp upstream of the 5′ region of ADRB1 and is in high LD with rs1801252, a missense variation in ADRB1 that codes for a change from serine to glycine at the 49th amino acid. According to HaploReg, rs12414657 (ADRB1) is associated with H3K4me3_Pro and 22_Promp methylation and acetylation histone modification in osteoblast primary cells, and this site is bound by Pol2, TAF1, or Pol24H8 in nonbone cell lines [70]. This suggests that this SNV is within the promoter region of ADRB1, as it is associated with H3K4me3, which is a histone modification that indicates a promoter region [97, 98]. This histone modification promotes chromatin remodeling that allows transcription factors to bind to that site [99], and Pol24H8 binding indicates a transcription factor binding site [100]. Thus, a variation at this location may affect the transcription of ADRB1. This is further supported in SCREEN, where rs12414657 (ADRB1) overlapped with a cis-regulatory element that expressed a cell-type agnostic proximal enhancer-like signature including high DNase, H3K4me4, H3K27ac, and CTCF markers (Z score > 1.64) [71]. These markers were not as strong in osteoblast primary cells where DNase-seq was not available, and H3K4me3 had the highest Z score at 1.08. Unlike rs12414657 (ADRB1), rs1801252 (ADRB1), a nonsynonymous SNV in high LD with rs12414657 (ADRB1), has been reported in association with resting heart rate and survival in patients with heart failure [101-108].

Sex-specific effects in β-adrenergic signaling have been shown previously. Specifically, β-adrenergic contractile response is greater in male mice than female mice [45], and female rabbit hearts have decreased capacity to respond to β-adrenergic stimulation as compared with male rabbit hearts [46]. In humans, sex differences related to β-adrenergic signaling are present with regard to susceptibility to heart failure, arrythmia, and other cardiovascular conditions, such as hypertension [45]. BBs are also reported to be less effective in women than men [49]. Osteoporosis and BMD also vary by sex, and osteoporosis is more prevalent in women, especially postmenopausal women, with ovariectomized mice serving as an in vivo model of postmenopausal osteoporosis [109]. Therefore, we performed sex-stratified genetic analyses and did not observe any overlap in findings between sexes. Sex-specific effects, if present, could reflect the differential occurrence of disease among sexes or could arise from differences in RNA expression, protein expression, or downstream response [42, 110-114]. The lack of replication of genetic findings across sexes may also be a result of power limitations, especially given that our findings in males did not replicate in external cohorts. Given that our top SNVs are noncoding SNVs, there may be a functional SNV in high LD with our top SNVs. Other possibilities include alteration in splicing efficiency, activation of cryptic splice sites, or altered expression of alternate transcripts. Furthermore, these SNVs may be involved in long-range gene regulation and influence the expression of remote genes as part of a regulatory element [115].

There are important limitations to our study. Our study is cross-sectional and thus our results are to be interpreted as an association, with causal mechanisms yet to be determined. Owing to limited power, we did not correct for multiple testing in the discovery cohort, but we did perform strict multiple testing correction in the validation cohorts using the Bonferroni method. We feel that these results should be validated in additional prospective studies for confirmation. In addition, because we did not perform a full GWAS because of limited power, there are many potential genes as well as long-range regions around our candidate genes that were not assessed that may have pharmacogenomics effects. Collider bias is also a potential concern as there are many genetic variations associated with height and BMI, which we used as covariates. However, we did not find that height or BMI had been previously associated with any of the SNVs submitted for validation or for SNVs in high LD with those SNVs. Additionally, for GCTA-COJO, Yang et al recommend a reference sample of more than 4000 individuals [69]. Although we did not have access to that sample size for our reference sample, previous simulated studies report an R2 greater than 0.9 when using a reference sample with more than 1000 individuals [69], which we did have. Another potential limitation of our findings is that the association of SNVs, miRNAs, and BB use could be due to confounding by treatment indication. We attempted to reduce the effects of confounding by conditioning on important clinical covariates and miRNA technical variables, but residual confounding is still possible. In addition, we performed a detailed series of sensitivity analyses for our top 2 SNVs in females in which we removed individuals being treated for bone disease or steroids, stratified by treatment with other medications known to have bone effects, examined the results of these SNVs at other skeletal sites, and looked at β1-selectivity and dose-specific models, and have found that the effect estimates were stable under all these scenarios. We were also underpowered to adjust for batch effects in our miRNA analysis, and our miRNA data, which are from whole blood, may not reflect expression in bone. However, these miRNAs have been previously associated with BMD, fractures, or osteoporosis using data from plasma, serum, or whole blood [36-38, 76-93], and we did test for potential confounding due to hemolysis. An additional limitation is the lack of diversity in the study sample, limiting our findings to participants of White, European ancestry. Also, these SNVs have not previously been cited in any GWAS study, which limits our knowledge to that obtained from our own data sets. We are also uncertain as to why these effects are different by sex, although prior evidence of sex differences in the effects of β-adrenergic signaling or power limitations may provide explanations. Another limitation is that while the study sample was homogenous, the validation cohorts used were heterogeneous in regards to sample design and demographics. This may have limited our power but strengthened the external validity of our findings.

To our knowledge, this is the first study to suggest an interaction between genes and BB use on BMD. We have identified 4 sex-specific genetic variations that map to or near genes involved in adrenergic signaling in bone and successfully validated 2 of them in external cohorts. Based on our previous work showing that BB use is associated with the presence of certain circulating miRNAs, we have also determined miRNAs associated with these SNVs and putative miRNA-mediated mechanisms by which these SNVs mediate the effect of BB use on BMD. We intend to validate these mechanisms in the future using in vivo, in vitro, and clinical models. In conclusion, our findings that BB associations with BMD may be modified by genetic variation suggest that studies evaluating the bone effects of BBs consider genetic variation in drug response.