Identifying Pleiotropic SNPs Associated With Femoral Neck and Heel Bone Mineral Density

Background: Genome-wide association studies (GWASs) routinely identify loci associated with risk factors for osteoporosis. However, GWASs with relatively small sample sizes still lack sufficient power to ascertain the majority of genetic variants with small to modest effect size, which may together truly influence the phenotype. The loci identified only account for a small percentage of the heritability of osteoporosis. This study aims to identify novel genetic loci associated with DXA-derived femoral neck (FNK) bone mineral density (BMD) and quantitative ultrasound of the heel calcaneus estimated BMD (eBMD), and to detect shared/causal variants for the two traits, to assess whether the SNPs or putative causal SNPs associated with eBMD were also associated with FNK-BMD.

Methods: Novel loci associated with eBMD and FNK-BMD were identified by the genetic pleiotropic conditional false discovery rate (cFDR) method. Shared putative causal variants between eBMD and FNK-BMD and putative causal SNPs for each trait were identified by the colocalization method. Mendelian randomization analysis addresses the causal relationship between eBMD/FNK-BMD and fracture.

Results: We identified 9,500 (cFDR < 9.8E-6), 137 (cFDR < 8.9E-4) and 124 SNPs associated with eBMD, FNK-BMD, and both eBMD and FNK-BMD, respectively, with 37 genomic regions where there was a SNP that influences both eBMD and FNK-BMD. Most genomic regions only contained putative causal SNPs associated with eBMD and 3 regions contained two distinct putative causal SNPs influenced both traits, respectively. We demonstrated a causal effect of FNK-BMD/eBMD on fracture.

Conclusion: Most of SNPs or putative causal SNPs associated with FNK-BMD were also associated with eBMD. However, most of SNPs or putative causal SNPs associated with eBMD were not associated with FNK-BMD. The novel variants we identified may help to account for the additional proportion of variance of each trait and advance our understanding of the genetic mechanisms underlying osteoporotic fracture.

Introduction

Osteoporosis is a complex disease that prevalently occurs in postmenopausal women and is diagnosed primarily by measuring bone mineral density (BMD) (Tella and Gallagher, 2014). Low femoral neck (FNK) BMD is a major risk factor of osteoporosis and remains the best predictor of primary osteoporotic fractures (Meng et al., 2018). It is established that there are differences of genetic determination among various skeletal sites (Yang et al., 2006). The heritability of BMD is estimated between 46 and 84% depending on the skeleton site studied (Mullin et al., 2018).

Although BMD can be measured in many ways, the “gold standard” diagnostic technique is dual energy absorptiometry (DXA) of the spine and hip (Tella and Gallagher, 2014). However, DXA is relatively expensive, and consequently, the genome-wide association studies (GWASs) of DXA-derived BMD include a relatively small size of samples, compromising the ability to detect genetic loci. Quantitative ultrasound (QUS) of the heel calcaneus estimated BMD (referred to as estimated BMD “eBMD” in this study) is quick, safe, relatively inexpensive and is cost-effective for GWASs with very large samples of individuals (Kemp et al., 2017). The correlations of eBMD with spine (r = 0.50) and hip BMD (r = 0.54) were modest but stronger in postmenopausal women (r = 0.60 for spine and r = 0.62 for hip, respectively) (Saadi et al., 2004). Furthermore, eBMD has been used in combination with other clinical risk factors to predict the risk of fracture (Khaw et al., 2004).

Genome-wide association studies is a powerful tool to identify susceptibility genetic variants of complex diseases. One previous GWAS identified 56 BMD-associated loci, accounting for approximately 5.8% of BMD variance (Estrada et al., 2012). The conditionally independent genome-wide significant lead single nucleotide polymorphisms (SNPs) account for 20% of the total variance of eBMD (Morris et al., 2019). However, despite the relative large sample sizes, GWASs with these sample sizes may still lack sufficient power to identify the majority of genetic variants with small to modest effect sizes, which may truly influence the phenotype (Pei et al., 2014).

Pleiotropy is a phenomenon in which a genetic variant affects multiple phenotypes. The presence of pleiotropy indicates that related traits may also share some genetic determinants (Greenbaum et al., 2017). Identifying potential shared causal variants among traits can help determine the overlapping etiologies of multifactorial disorders. Therefore, we considered variants related to both BMD and normal bone physiology. There have been limited efforts to ascertain the mechanisms that alter bone physiology or the causal genes underlying the GWAS loci (Sabik and Farber, 2017). GWASs allow us to obtain the estimates of effect size for all those genetic variants, so it is possible to estimate shared genetic determinants by checking the correlations between the effect sizes across traits/bone parameters, which does not need to measure multiple traits in the same individuals (Bulik-Sullivan et al., 2015). We used GWAS-derived summary statistics to identify the shared genetic variants for two related phenotypes with the genetic pleiotropy-informed conditional false discovery rate (cFDR) method (Andreassen et al., 2013). This method incorporates the GWAS-derived summary statistics for two related phenotypes to test the association of genetic variants with one phenotype conditional on the different strengths of the association of these variants with the second phenotype (Greenbaum et al., 2017). cFDR improves the statistical power and identifies additional variants associated with each trait and pleiotropic variants from the studies of traits or diseases (Lv et al., 2017). Bayesian colocalization models identify genomic regions containing the same association signal by GWAS and colocalizing expression quantitative trait loci (eQTL) signals or GWAS signals of two traits (Giambartolomei et al., 2014; Pickrell et al., 2016). This model is the basis to integrally consider the effect sizes of all genetic variants in the pre-defined regions (e.g., approximately independent linkage disequilibrium (LD) blocks) (Berisa and Pickrell, 2016) and calculate the posterior probabilities of each region containing one causal variant associated with both traits (Giambartolomei et al., 2014; Pickrell et al., 2016).

In this study, we aim to compare those genetic loci associated with DXA-derived FNK-BMD or eBMD and putative causal SNPs associated with each trait. Novel variants with pleiotropic or causal effects on eBMD and FNK-BMD were identified by cFDR and colocalization analysis of GWAS summary statistics. Moreover, we investigated the relationship between fracture and FNK-BMD/eBMD by Mendelian randomization (MR) analysis of summary statistics in order to assess if FNK-BMD/eBMD causally underlie fracture risk.

Materials and Methods

GWAS Datasets

The summary statistics analyzed were derived from three publicly available online GWAS datasets: FNK-BMD dataset, eBMD dataset and fracture dataset. The FNK-BMD GWAS dataset contained the association results for approximately 10 million SNPs and is derived from a GWAS meta-analysis of several studies with 32,965 European subjects published by the Genetic Factors for Osteoporosis (GEFOS) Consortium (Zheng et al., 2015). A GWAS meta-analysis was conducted in 426,824 subjects of European ancestry to study the associations of over 14 million SNPs with eBMD (Morris et al., 2019). The fracture dataset was generated from a GWAS meta-analysis of 264,973 European subjects, containing summary statistics of more than 13 million SNPs (Morris et al., 2019).

Data Preparation

We pruned SNPs in the FNK-BMD and eBMD datasets before cFDR analysis using the genotype data from the 1000 Genomes project as a reference. Briefly, two datasets were combined to select overlapping SNPs, and then pairs of SNPs with a high correlation were pruned with a LD-based algorithm (Greenbaum et al., 2017). Default values of the PLINK 1.9 software (50, 5, 0.2) were set as parameters when calculating LD values (r²) between each SNP pair (Greenbaum et al., 2017). The LD was computed for windows which contained 50 SNPs. The SNP with lower frequency of the minor allele was excluded for each pair with r² > 0.2. The calculation window was then shifted forward by 5 SNPs. Then repeated above process until each pair of SNPs were in low LD. Finally, 1,783,638 SNPs were included in the following analyses.

Data Analysis

Estimating Heritability and Genetic Correlation

Because the previous FNK-BMD GWAS did not show the proportion of variance explained by the significant findings in their study, the heritability explained by the previous GWAS results was estimated using the GWAS summary statistics. The stratified LD score regression partitioned heritability and estimated the genetic correlation between FNK-BMD and eBMD (Finucane et al., 2015).

Conditional FDR Calculation and Conditional Manhattan Plot

Conditional false discovery rate is an extension of the FDR method, which combines the GWAS summary statistics of two phenotypes to evaluate the probability of the association of genetic variants with the principal phenotype conditioned on the strength of the association with the conditional phenotype (Liu et al., 2018). Detailed principles and formulas of the method were described earlier by Greenbaum (Greenbaum et al., 2017). P-values from GWAS-derived summary statistics of eBMD and FNK-BMD were searched to identify common SNPs, through SNP pruning for cFDR and conjunction cFDR (ccFDR) analyses. Then we calculated the cFDR for each SNP where FNK-BMD is the principal phenotype conditioned on the strength of association with eBMD (FNK-BMD | eBMD) and vice versa (eBMD | FNK-BMD).

Additionally, ccFDR value was computed to identify genetic variants associated with both eBMD and FNK-BMD simultaneously. The ccFDR value is defined as the probability that the given SNP was wrongly identified to be associated with both phenotypes. The thresholds of cFDR were set at the maximum cFDR with p value less than 5 × 10^–8 in the principal phenotype (Liley and Wallace, 2015). Manhattan plots were created with R software.

Pleiotropic Enrichment Estimation

To assess pleiotropy between eBMD and FNK-BMD, conditional QQ and fold-enrichment were plotted. We presented the QQ curves of the quantiles of nominal –log₁₀ (p)-values obtained from GWAS summary statistics for the association of the SNPs subsets in the conditional phenotype below each significance threshold. The x-axis and y-axis are empirical quantiles of the nominal p values (empirical distribution functions) and nominal p values (-log₁₀ (p)), respectively. Herein, nominal -log₁₀ (p)-values were stratified by p-values of conditional phenotype with the cutoffs as p < 1 (expected base line, all SNPs), p < 0.1, p < 0.01, p < 0.001 and p < 0.0001. A leftward curve defection from the expected base line in the QQ plots represents pleiotropic enrichment shared by two phenotypes.

Nominal p values (- log₁₀ (p)) and fold enrichments were plotted on the x-axis and the y-axis in fold-enrichment plots, respectively. We presented fold-enrichment plots of nominal - log₁₀ (p)-values for eBMD SNPs below the standard GWAS threshold of p < 5 × 10^–8 and for subsets of SNPs determined by the significance of their association with FNK-BMD and vice versa. A shift above the expected baseline indicates pleiotropy between the two phenotypes (Schork et al., 2013).

Colocalization Analysis

The Bayesian colocalization method estimates the probability of each genomic locus containing a genetic variant that is causally associated with both traits (Guo et al., 2015). The Bayesian approach summarizes the evidence of each genetic region for five mutually exclusive hypotheses simultaneously:

H₀: There exists no SNPs in the region that are associated with either trait,

H₁ (Model 1): There exists one causal SNP in the region that is associated with the first trait,

H₂ (Model 2): There exists a causal SNP in the region that is associated with the second trait,

H₃ (Model 3): There exists a single causal SNP in the region that are associated with both traits,

H₄ (Model 4): There exist two distinct causal SNPs in the region, one for each trait.

Support for each hypothesis is quantified on the basis of posterior probabilities (PP), defined by PP₀, PP₁, PP₂, PP₃, or PP₄ (Guo et al., 2015).

Herein, colocalization was analyzed with gwas-pw software (Pickrell et al., 2016), which is available at https://github.com/joepickrell/gwas-pw. The method requires the effect size estimates and standard errors of each SNP for the two phenotypes. The whole genome was split into 1,703 non-overlapping regions (blocks). The regional posterior probability of the region containing one SNP causally associated with both phenotypes larger than 0.9 corresponded with a FDR < 0.1. The SNP in each region with the largest posterior probability was defined as putative causal SNP.

Genes and SNPs Annotation

ANNOVAR software was employed to annotate the SNPs of interest (Alexander et al., 2017). For genes of interest, Gene Ontologies (GO) and Kyoto Encyclopedia of Genes and Genomes (KEEG) pathways were analyzed using DAVID 6.7¹.

Results

Heritability and Genetic Correlation

The previous FNK-BMD GWAS accounted for 9.1% of the SNP-based heritability (95% CI: 8.1–10.1%). The genetic correlation between FNK-BMD and eBMD was estimated to be 0.64 (p = 5.9 × 10^–61).

Pleiotropy Between eBMD and FNK-BMD

The conditional QQ plot for eBMD conditional on FNK-BMD indicated enrichment of pleiotropic SNPs across the different significance thresholds for FNK-BMD (Figure 1). An obvious upward shift from the expected baseline was detected when restricting the SNPs subgroup with a stronger level of association with the conditional phenotype, indicating an increase in the number of true associations for a given FNK-BMD p-value (Figure 1A). Similar enrichment of pleiotropic SNPs for FNK-BMD given eBMD, was evidenced by a similar departure pattern across the different curves (Figure 1B). For any given eBMD nominal p-value, a deviation from the null line demonstrated a greater proportion of true pleiotropic associations.

FIGURE 1

Figure 1. Stratified QQ (A,B) and Enrichment (C,D) plots. Stratified QQ plots of nominal versus empirical –log₁₀ p-values in eBMD as a function of significance of the association with FNK-BMD (A), and in FNK-BMD as a function of significance of the association with eBMD (B). Fold-enrichment plots (C,D) of enrichment versus nominal –log₁₀ p-values for eBMD below the standard GWAS threshold of p < 5 × 10^–8 as a function of significance of the association with FNK-BMD, and FNK-BMD below the standard GWAS threshold of p < 5 × 10^–8 as a function of significance of the association with eBMD.

The fold-enrichment plots showed SNP enrichment for eBMD with FNK-BMD across different significant levels, and also for FNK-BMD with eBMD (Figures 1C,D). It illustrated that between the SNPs subgroup with the strongest association in the conditional phenotype and all SNPs groups, the proportion of SNPs associated with the principal phenotype increased by approximately 30-fold at the genome-wide significance level.

eBMD and FNK-BMD Loci Detected by cFDR

We identified 9,500 SNPs significantly associated with eBMD (cFDR < 9.8E-06), given their associations with FNK-BMD (Supplementary Table S1). All of these SNPs had p values less than 1 × 10^–5 and 8,933 reached the genome-wide significance at p < 5 × 10^–8 in the original GWAS meta-analysis for eBMD (Kemp et al., 2017). Among the 9,500 SNPs we identified, 50.2% (4,769 SNPs) were intronic, 1.7% (n = 163) were juxtaposed to the 3’ or 5’ terminals of corresponding genes, and 0.7% (n = 71) were in the 3’ or 5’ UTR (Supplementary Table S1).

Conditional on the association of SNPs with eBMD, we identified 137 significant SNPs (cFDR < 8.9E-04) associated with FNK-BMD (Supplementary Table S2). In the original FNK-BMD GWAS meta-analysis (Zheng et al., 2015), 7 SNPs had p values less than 1 × 10^–5 while 4 of them reached the standard genome-wide significance. Among the 137 SNPs identified, 56.2% (77 SNPs) were intronic, and 1.6% were juxtaposed to the 3’ or 5’ terminals of corresponding genes (Supplementary Table S2).

Pleiotropic SNPs for Both eBMD and FNK-BMD

We identified 124 pleiotropic SNPs that were significantly associated with both eBMD and FNK-BMD using ccFDR method (ccFDR < 8.9E-04) (Supplementary Table S3 and Figure 2). It suggested that most of SNPs (90.5%, $\frac{124}{137}$) associated with FNK-BMD were also associated with eBMD. However, most of SNPs (98.7%, 1-$\frac{124}{9500}$) associated with eBMD were not associated with FNK-BMD. Although, this can be partly explained by the relatively large sample size of eBMD study and the relatively small sample size of FNK-BMD study. Among these SNPs, 25 previously reached the standard GWAS significance (p < 5 × 10^–8) both in the original eBMD and FNK-BMD GWASs. Another 95 pleiotropic SNPs were significantly associated with eBMD in original eBMD GWAS (original p < 5 × 10^–8). Four (rs79676715, rs1381635, rs7816021 and rs4263799) of these statistically significant SNPs identified by ccFDR were not reported by the original eBMD and FNK-BMD GWASs.

FIGURE 2

Figure 2. Conjunction Manhattan plot of conjunction –log₁₀ FDR values for eBMD and FNK-BMD. The red line marking the conditional –log₁₀ ccFDR value of 3.1 corresponds to a conjunction FDR < 8.9E-04.

The 124 pleiotropic SNPs were located within 73 different genes dispersed over 18 autosomes. GO analyses indicated that these genes were significantly enriched in the canonical Wnt signaling pathway (GO: 0060070, p = 3.8E-06) and positive regulation of transcription, DNA-templated (GO: 0045893, p = 1.3E-05) (Table 1). KEGG pathway analysis confirmed that pleiotropic genes were significantly enriched in Wnt signaling pathway (hsa04310, p = 1.8E-06), and also identified signaling pathways regulating pluripotency of stem cells (hsa04550, p = 3.9E-05) and genes associated with basal cell carcinoma (hsa05217, p = 1.5E-05) (Table 2).

TABLE 1

Table 1. GO terms enriched by pleiotropic genes.

TABLE 2

Table 2. KEGG pathways enriched by pleiotropic genes.

Colocalization Analysis

We identified 37 genomic regions containing a SNP that influences both eBMD and FNK-BMD (model 3) (Supplementary Table S4), 3 genomic regions containing two distinct causal SNPs associated with both traits, respectively (model 4) and 627 regions containing causal SNPs only for eBMD (model 2) (Supplementary Table S5), at a threshold of the posterior probability of model 3 larger than 0.9. The 37 putative causal SNPs annotated to 53 different genes, which were significantly enriched in the biological processes of canonical Wnt signaling (GO: 0060070, p = 6.5E-07), positive regulation of transcription (GO: 0045893, p = 6.7E-06), and DNA-templated and positive regulation of transcription from RNA polymerase II promoter (GO: 0045944, p = 1.6E-05) (Table 3). Block 42 (Genomic region: chr1: 66939404–68477895), containing the candidate gene GNG12-AS1(WLS), was associated with both eBMD and FNK-BMD, and the putative casual SNP (rs2566759) in this region is associated with spine BMD in a previous study (Mullin et al., 2016). Of the 37 genomic regions identified, 26 contained pleiotropic SNPs that were identified in ccFDR analysis. We noted that the putative causal SNPs identified by colocalization analysis were distinct from pleiotropic SNPs in each genomic region, suggesting the ccFDR method is apt to identify the variant associated with both phenotypes instead of the causal one.

TABLE 3

Table 3. GO terms enriched by causal genes.

Discussion

In the present study, we combined two GWAS summary statistics by cFDR to identify the pleiotropic and/or novel SNPs associated with eBMD and/or FNK-BMD. We discovered 9,500 SNPs and 137 SNPs for eBMD and FNK-BMD, respectively. It suggested that the majority of SNPs associated with eBMD were not identified for FNK-BMD. Conjunction cFDR analysis identified 124 pleiotropic SNPs that were significantly associated with both eBMD and FNK-BMD. We also identified 37 putative causal SNPs that may underlie both phenotypes using colocalization analysis. It suggested that both BMD phenotype measurements shared some causal genetic variants in some regions, but in most of independent regions (94.5%, $\frac{3+627}{37+3+627}$) the putative causal SNPs for both phenotypes were different.

As osteoporosis is a complex disease with a strong genetic component, scientists have strived to identify SNPs associated with BMD and/or fracture risk. In the present study, we identified 37 candidate SNPs which are causally associated with both eBMD and FNK-BMD. Among these SNPs, SNP rs884205 is located in the 3’UTR of tumor necrosis factor (TNF) receptor super family member 11 a (TNFRSF11A), a gene that is a key driver for BMD and fracture risk (Tu et al., 2015). TNFRSF11A encodes RANK, a receptor that regulates bone remodeling. Genetic variation within the TNFRSF11A gene is associated with BMD in several previous GWASs (Styrkarsdottir et al., 2009). This SNP is associated with osteoporotic fractures in a Chinese population (Guo et al., 2012).

Dozens of studies in pre- and postmenopausal women, as well as in men, reported numerous of genes that are responsible for regulating serum periostin levels and bone microstructure (Pepe et al., 2018). These include genes with polymorphisms in low-density lipoprotein receptor-related protein (LRP5), estrogen-receptor 1 (ESR1) and tumor necrosis factor receptor super family member 11 (TNFRSF11) (Pepe et al., 2018). Several novel putative causal SNPs were near those genes identified in a previous study, including rs11228240 (LRP5), rs2941741 (ESR1), and rs78667121 (TNFRSF11).

The SNP rs2566752 is an intronic variant located in the wntless Wnt ligand secretion mediator (WLS) gene, a known BMD-associated locus that encodes an integral component of the Wnt ligand secretion pathway (Mullin et al., 2016). The Wnt proteins are members of an evolutionarily conserved family of secreted signaling molecules and are critical for prenatal and postnatal bone development (de Groot et al., 2014; Zhong et al., 2015). Variants within the WLS gene are associated with BMD in several previous GWASs (Rivadeneira et al., 2009). A large meta-analysis also identified rs2566752 as the maximally associated variant from the WLS gene region for both lumbar spine and FNK-BMD (Mullin et al., 2016). Interestingly, bone geometric parameters and BMD pleiotropy were previously detected at the WLS locus (Roshandel et al., 2011). These findings suggest that the associations between rs2566752 and BMD are likely mediated through regulatory effects on the WLS gene.

Among the 37 putative causal SNPs, five SNPs overlap with enhancer- or promoter-type chromatin in osteoblasts (i.e., rs9594738, rs7587430, rs11023993, rs4776341, and rs71390846) (Kundaje et al., 2015). Interestingly, rs9594738, overlaps a weak enhancer-type chromatin, lies ∼54.7 kb downstream of A-kinase anchor protein 11 (AKAP11) and ∼184.7 kb upstream of TNFSF11, or commonly known as RANKL. The role of APAK11 in bone biology is unclear, and it has broad expression in many types of cell and tissue (Encode Project Consortium, 2012). However, RANKL is well known to be involved in RANKL/RANK/osteoprotegerin signaling, which is one of the key signaling pathways controlling bone resorption and formation (Chen et al., 2018). In a recent study, the CRISPR/Cas9 deletion of rs9594738 and a pleiotropic SNP rs9533090 was shown to reduce RANKL expression at both mRNA level and protein level in U2-OS cells (Zhu et al., 2018).

Another putative causal SNP, rs4776341, overlaps with an active enhancer-type chromatin in osteoblasts, is an intronic variant located in the sterile alpha motif domain containing 3 (SMAD3) gene, which is part of the important TGF-β1 signaling pathway for osteoblast and osteoclast differentiation via the TGF-β1 signaling pathway (Yang et al., 2017). Interestingly, the SNP rs4776341 is associated with eBMD and FNK-BMD by the cFDR method. However, in most causal regions, the putative causal SNPs identified by the colocalization analysis were different from the pleiotropic SNPs identified by the cFDR method. This may be because the real causal SNP was removed using an LD-based pruning method, which was employed to remove one SNP of pairs with R² value greater than 0.2 during the process of data preparation in the cFDR method. For example, the pleiotropic SNP rs6426749 and putative causal SNP rs34553872 in block 14 (chr1:21736983-23086441) were in strong LD (r² = 1.0). The putative causal SNP rs34553872 was removed during SNP pruning.

Investigating a genetic variant influencing multiple phenotypes is informative in several applications. It can uncover the shared genetic mechanisms between closely related phenotypes, as well as the molecular function of a gene (Solovieff et al., 2013). The cFDR approach was applied herein to identify pleiotropic genetic variants. This method is based on the concept that a SNP with a significant effect on one trait is more likely to be associated with another correlated trait, and therefore is more likely to be detected by combining multiple independent studies (Greenbaum et al., 2017). This approach can not only improve the statistical power for detecting novel variants with true associations, but can detect small to intermediate effect sizes (Andreassen et al., 2013). The genetic enrichment observed in conditional QQ plots illustrates that the power of detecting true association effects greatly enhanced as the cFDR value for a given nominal p value decreased. Andreassen and colleagues showed that the cFDR method resulted 15–20-fold increase in the number of SNPs under the same FDR threshold of 0.05, compared with the traditional unconditional FDR method (Andreassen et al., 2013). The cFDR method can identify genetic variant associated with the “principal” phenotype, as long as the p-value of association with both the “principal” and the “conditional” phenotypes is below a certain p-value threshold (Andreassen et al., 2013; Hackinger and Zeggini, 2017). This method can be used as a variant prioritization tool to identify other variants associated with one phenotype by utilizing information of the association with the second phenotype (Hackinger and Zeggini, 2017). In addition, the cross-phenotype effect can be detected by calculating the maximum value of cFDR for two traits.

The colocalization method is based on Z scores of the SNPs from two input studies, testing an excess of shared signals, and identifying a set of variants with evidence for the association with the two traits (Hackinger and Zeggini, 2017). Additionally, we investigated whether causal SNPs are shared by combining information of multiple SNPs in a region by a colocalization analysis (Fortune et al., 2015). The main strength of this method lies in its speed and analytical forms, and the fact that it can use variants’ effect size estimates and standard errors when they are the only data available (Giambartolomei et al., 2014; Pickrell et al., 2016). This approach focuses on the effects of all SNPs, including those that do not reach genome-wide significance (Fortune et al., 2015). Thus, variants of small effect sizes are more readily identified.

The 124 pleiotropic SNPs identified by the ccFDR analysis were located in 37 independent blocks which were defined in colocalization analysis. Among these blocks, 26 overlapped with causal regions identified by colocalization analysis. This analysis implies that the cFDR method is more powerful for detecting pleiotropic SNPs. Although the purpose of colocalization analysis was to identify the putative causal SNPs affecting both phenotypes, most of the pleiotropic regions (belonging to model 3 and model 4) identified by colocalization analysis contained a putative causal SNP affecting both phenotypes, rather than two distinct putative causal SNPs in the region (one for each trait). Both the cFDR and colocalization results suggest that there is a significant genetic correlation between the two phenotypes. Consistent with the results of the genetic correlation analyses, SNPs influencing FNK-BMD were strongly correlated with eBMD (genetic correlation 0.64, p = 5.92E-61). The genetic correlation analyses between two traits indicated that the majority of shared genetic determinants remains to be discovered. The causal association between BMD and fracture risk is well established. Herein, we also found that both FNK-BMD and eBMD have a causal effect on fracture.

Estimated BMD is derived based on the equation: 0.002592 × (BUA+SOS)-3.687. Previously eBMD study indicated that eBMD is not an actual measurement of calcaneal BMD (Frost et al., 2000). Although eBMD is highly correlated with DXA-derived BMD, there is still some SNPs which is associated with DXA-derived BMD but not associated with eBMD. Therefore, we concluded that eBMD can be a good proxy of BMD but cannot be a replacement of DXA-derived BMD, especially used for diagnosis.

We suspected the differences in GWAS sample size could have had an impact on the cFDR results. Previous study has performed a large scale GWAS on total body BMD (TB-BMD) and showed that TB-BMD could be a good proxy of FNK-BMD. We performed cFDR method to identify novel variants with pleiotropic effects on eBMD and TB-BMD using TB-BMD GWAS summary statistics which comprised over 66,000 individuals (Medina-Gomez et al., 2018). The results showed that most of SNPs associated with eBMD were not associated with TB-BMD. The cFDR results are mainly affected by the rank of each SNP’s p value, not the exact p value in the respective GWAS studies of the two traits. Although this can be partly explained by the different sample sizes of the TB-BMD and eBMD datasets, the results should be reasonable and objective when the studied GWASs datasets are used.

In summary, by applying cFDR and colocalization methods, we identified pleiotropic SNPs associated with two risk traits (e.g., eBMD and FNK-BMD) of osteoporotic fracture, and putative causal SNPs that influence both traits. These findings suggested that eBMD cannot be used for a replacement of DXA-derived BMD (e.g., FNK-BMD). These findings may offer novel pathophysiological insight, and uncover potential targets and pathways for the trait-associated loci impacting bone health in osteoporosis.

Data Availability Statement

The summary statistics analyzed were derived from three publicly available online GWAS datasets: FNK-BMD dataset, eBMD dataset, and fracture dataset.

Author Contributions

H-WD, PH, and X-HM: study design. PH, X-HM, and XL: data collection and analysis. PH and XZ: drafting manuscript. H-WD, MS, X-HM, QZ, F-YD, and R-LJ: revising manuscript content. PH and X-HM took responsibility for the integrity of the data analysis. All authors contributed to the preparation and approval of the final manuscript.

Funding

H-WD was partially supported by grants from NIH (R01 AR069055, U19 AG055373, R01 MH104680, R01 AR059781, and P20 GM109036), and the Edward G. Schlieder Endowment. MS was supported by a grant from NIH P20 GM121325 and a Prabhu Endowed Professorship.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

We thank all the study subjects for volunteering to participate in the study. We thank Genetic Factors for Osteoporosis Consortium provided the BMD dataset to download.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fgene.2020.00772/full#supplementary-material

TABLE S1 | eBMD loci identified with cFDR.

TABLE S2 | FNK-BMD loci identified with cFDR.

TABLE S3 | SNPs identified by ccFDR analysis.

TABLE S4 | Colocalization analysis identified 37 significant genomic regions.

TABLE S5 | Results for Model 2 and 4.

Footnotes

1. ^ https://david.ncifcrf.gov/summary.jsp

References

Alexander, J., Mantzaris, D., Georgitsi, M., Drineas, P., and Paschou, P. (2017). Variant Ranker: a web-tool to rank genomic data according to functional significance. BMC Bioinformatics 18:341. doi: 10.1186/s12859-017-1752-3

PubMed Abstract | CrossRef Full Text | Google Scholar

Andreassen, O. A., Thompson, W. K., Schork, A. J., Ripke, S., Mattingsdal, M., Kelsoe, J. R., et al. (2013). Improved detection of common variants associated with schizophrenia and bipolar disorder using pleiotropy-informed conditional false discovery rate. PLoS Genet. 9:e1003455. doi: 10.1371/journal.pgen.1003455

PubMed Abstract | CrossRef Full Text | Google Scholar

Berisa, T., and Pickrell, J. K. (2016). Approximately independent linkage disequilibrium blocks in human populations. Bioinformatics 32, 283–285.

Google Scholar

Bulik-Sullivan, B., Finucane, H. K., Anttila, V., Gusev, A., Day, F. R., Loh, P. R., et al. (2015). An atlas of genetic correlations across human diseases and traits. Nat. Genet. 47, 1236–1241. doi: 10.1038/ng.3406

PubMed Abstract | CrossRef Full Text | Google Scholar

Chen, X., Wang, Z., Duan, N., Zhu, G., Schwarz, E. M., and Xie, C. (2018). Osteoblast-osteoclast interactions. Connect Tissue Res. 59, 99–107.

Google Scholar

de Groot, R. E., Rappel, S. B., Lorenowicz, M. J., and Korswagen, H. C. (2014). Protein kinase CK2 is required for Wntless internalization and Wnt secretion. Cell Signal. 26, 2601–2605. doi: 10.1016/j.cellsig.2014.08.016

PubMed Abstract | CrossRef Full Text | Google Scholar

Encode Project Consortium (2012). An integrated encyclopedia of DNA elements in the human genome. Nature 489, 57–74. doi: 10.1038/nature11247

PubMed Abstract | CrossRef Full Text | Google Scholar

Estrada, K., Styrkarsdottir, U., Evangelou, E., Hsu, Y. H., Duncan, E. L., Ntzani, E. E., et al. (2012). Genome-wide meta-analysis identifies 56 bone mineral density loci and reveals 14 loci associated with risk of fracture. Nat. Genet. 44, 491–501.

Google Scholar

Finucane, H. K., Bulik-Sullivan, B., Gusev, A., Trynka, G., Reshef, Y., Loh, P. R., et al. (2015). Partitioning heritability by functional annotation using genome-wide association summary statistics. Nat. Genet. 47, 1228–1235. doi: 10.1038/ng.3404

PubMed Abstract | CrossRef Full Text | Google Scholar

Fortune, M. D., Guo, H., Burren, O., Schofield, E., Walker, N. M., Ban, M., et al. (2015). Statistical colocalization of genetic risk variants for related autoimmune diseases in the context of common controls. Nat. Genet. 47, 839–846. doi: 10.1038/ng.3330

PubMed Abstract | CrossRef Full Text | Google Scholar

Frost, M. L., Blake, G. M., and Fogelman, I. (2000). Does quantitative ultrasound imaging enhance precision and discrimination? Osteoporos. Int. 11, 425–433. doi: 10.1007/s001980070110

PubMed Abstract | CrossRef Full Text | Google Scholar

Giambartolomei, C., Vukcevic, D., Schadt, E. E., Franke, L., Hingorani, A. D., Wallace, C., et al. (2014). Bayesian test for colocalisation between pairs of genetic association studies using summary statistics. PLoS Genet. 10:e1004383. doi: 10.1371/journal.pgen.1004383

PubMed Abstract | CrossRef Full Text | Google Scholar

Greenbaum, J., Wu, K., Zhang, L., Shen, H., Zhang, J., and Deng, H. W. (2017). Increased detection of genetic loci associated with risk predictors of osteoporotic fracture using a pleiotropic cFDR method. Bone 99, 62–68. doi: 10.1016/j.bone.2017.03.052

PubMed Abstract | CrossRef Full Text | Google Scholar

Guo, H., Fortune, M. D., Burren, O. S., Schofield, E., Todd, J. A., and Wallace, C. (2015). Integration of disease association and eQTL data using a Bayesian colocalisation approach highlights six candidate causal genes in immune-mediated diseases. Hum. Mol. Genet. 24, 3305–3313. doi: 10.1093/hmg/ddv077

PubMed Abstract | CrossRef Full Text | Google Scholar

Guo, Y., Wang, J. T., Liu, H., Li, M., Yang, T. L., Zhang, X. W., et al. (2012). Are bone mineral density loci associated with hip osteoporotic fractures? A validation study on previously reported genome-wide association loci in a Chinese population. Genet. Mol. Res. 11, 202–210. doi: 10.4238/2012.january.31.1

PubMed Abstract | CrossRef Full Text | Google Scholar

Hackinger, S., and Zeggini, E. (2017). Statistical methods to detect pleiotropy in human complex traits. Open Biol. 7:170125. doi: 10.1098/rsob.170125

PubMed Abstract | CrossRef Full Text | Google Scholar

Kemp, J. P., Morris, J. A., Medina-Gomez, C., Forgetta, V., Warrington, N. M., Youlten, S. E., et al. (2017). Identification of 153 new loci associated with heel bone mineral density and functional involvement of GPC6 in osteoporosis. Nat. Genet. 49, 1468–1475. doi: 10.1038/ng.3949

PubMed Abstract | CrossRef Full Text | Google Scholar

Khaw, K. T., Reeve, J., Luben, R., Bingham, S., Welch, A., Wareham, N., et al. (2004). Prediction of total and hip fracture risk in men and women by quantitative ultrasound of the calcaneus: EPIC-Norfolk prospective population study. Lancet 363, 197–202. doi: 10.1016/s0140-6736(03)15325-1

CrossRef Full Text | Google Scholar

Kundaje, A., Meuleman, W., Ernst, J., Bilenky, M., Yen, A., Heravi-Moussavi, A., et al. (2015). Integrative analysis of 111 reference human epigenomes. Nature 518, 317–330.

Google Scholar

Liley, J., and Wallace, C. (2015). A pleiotropy-informed Bayesian false discovery rate adapted to a shared control design finds new disease associations from GWAS summary statistics. PLoS Genet. 11:e1004926. doi: 10.1371/journal.pgen.1004926

PubMed Abstract | CrossRef Full Text | Google Scholar

Liu, H. M., He, J. Y., Zhang, Q., Lv, W. Q., Xia, X., Sun, C. Q., et al. (2018). Improved detection of genetic loci in estimated glomerular filtration rate and type 2 diabetes using a pleiotropic cFDR method. Mol. Genet. Genomics 293, 225–235. doi: 10.1007/s00438-017-1381-6

PubMed Abstract | CrossRef Full Text | Google Scholar

Lv, W. Q., Zhang, X., Zhang, Q., He, J. Y., Liu, H. M., Xia, X., et al. (2017). Novel common variants associated with body mass index and coronary artery disease detected using a pleiotropic cFDR method. J. Mol. Cell. Cardiol. 112, 1–7. doi: 10.1016/j.yjmcc.2017.08.011

PubMed Abstract | CrossRef Full Text | Google Scholar

Medina-Gomez, C., Kemp, J. P., Trajanoska, K., Luan, J., Chesi, A., Ahluwalia, T. S., et al. (2018). Life-course genome-wide association study meta-analysis of total body BMD and assessment of age-specific effects. Am. J. Hum. Genet. 102, 88–102.

Google Scholar

Meng, X. H., Shen, H., Chen, X. D., Xiao, H. M., and Deng, H. W. (2018). Inferring causal relationships between phenotypes using summary statistics from genome-wide association studies. Hum. Genet. 137, 247–255. doi: 10.1007/s00439-018-1876-1

PubMed Abstract | CrossRef Full Text | Google Scholar

Morris, J. A., Kemp, J. P., Youlten, S. E., Laurent, L., Logan, J. G., Chai, R. C., et al. (2019). An atlas of genetic influences on osteoporosis in humans and mice. Nat. Genet. 51, 258–266.

Google Scholar

Mullin, B. H., Walsh, J. P., Zheng, H. F., Brown, S. J., Surdulescu, G. L., Curtis, C., et al. (2016). Genome-wide association study using family-based cohorts identifies the WLS and CCDC170/ESR1 loci as associated with bone mineral density. BMC Genomics 17:136. doi: 10.1186/s12864-016-2481-0

PubMed Abstract | CrossRef Full Text | Google Scholar

Mullin, B. H., Zhu, K., Xu, J., Brown, S. J., Mullin, S., Tickner, J., et al. (2018). Expression quantitative trait locus study of bone mineral density GWAS variants in human osteoclasts. J. Bone Miner. Res. 33, 1044–1051. doi: 10.1002/jbmr.3412

PubMed Abstract | CrossRef Full Text | Google Scholar

Pei, Y. F., Zhang, L., Papasian, C. J., Wang, Y. P., and Deng, H. W. (2014). On individual genome-wide association studies and their meta-analysis. Hum. Genet. 133, 265–279. doi: 10.1007/s00439-013-1366-4

PubMed Abstract | CrossRef Full Text | Google Scholar

Pepe, J., Bonnet, N., Herrmann, F. R., Biver, E., Rizzoli, R., Chevalley, T., et al. (2018). Interaction between LRP5 and periostin gene polymorphisms on serum periostin levels and cortical bone microstructure. Osteoporos. Int. 29, 339–346. doi: 10.1007/s00198-017-4272-0

PubMed Abstract | CrossRef Full Text | Google Scholar

Pickrell, J. K., Berisa, T., Liu, J. Z., Segurel, L., Tung, J. Y., and Hinds, D. A. (2016). Detection and interpretation of shared genetic influences on 42 human traits. Nat. Genet. 48, 709–717. doi: 10.1038/ng.3570

PubMed Abstract | CrossRef Full Text | Google Scholar

Rivadeneira, F., Styrkarsdottir, U., Estrada, K., Halldorsson, B. V., Hsu, Y. H., Richards, J. B., et al. (2009). Twenty bone-mineral-density loci identified by large-scale meta-analysis of genome-wide association studies. Nat. Genet. 41, 1199–1206. doi: 10.1038/ng.446

PubMed Abstract | CrossRef Full Text | Google Scholar

Roshandel, D., Thomson, W., Pye, S. R., Boonen, S., Borghs, H., Vanderschueren, D., et al. (2011). Polymorphisms in genes involved in the NF-kappaB signalling pathway are associated with bone mineral density, geometry and turnover in men. PLoS One 6:e28031. doi: 10.1371/journal.pone.0028031

PubMed Abstract | CrossRef Full Text | Google Scholar

Saadi, H. F., Reed, R. L., Carter, A. O., and Al-Suhaili, A. R. (2004). Correlation of quantitative ultrasound parameters of the calcaneus with bone density of the spine and hip in women with prevalent hypovitaminosis D. J. Clin. Densitom. 7, 313–318. doi: 10.1385/jcd:7:3:313

CrossRef Full Text | Google Scholar

Sabik, O. L., and Farber, C. R. (2017). Using GWAS to identify novel therapeutic targets for osteoporosis. Transl. Res. 181, 15–26. doi: 10.1016/j.trsl.2016.10.009

PubMed Abstract | CrossRef Full Text | Google Scholar

Schork, A. J., Thompson, W. K., Pham, P., Torkamani, A., Roddey, J. C., Sullivan, P. F., et al. (2013). All SNPs are not created equal: genome-wide association studies reveal a consistent pattern of enrichment among functionally annotated SNPs. PLoS Genet. 9:e1003449. doi: 10.1371/journal.pgen.1003449

PubMed Abstract | CrossRef Full Text | Google Scholar

Solovieff, N., Cotsapas, C., Lee, P. H., Purcell, S. M., and Smoller, J. W. (2013). Pleiotropy in complex traits: challenges and strategies. Nat. Rev. Genet. 14, 483–495. doi: 10.1038/nrg3461

PubMed Abstract | CrossRef Full Text | Google Scholar

Styrkarsdottir, U., Halldorsson, B. V., Gretarsdottir, S., Gudbjartsson, D. F., Walters, G. B., Ingvarsson, T., et al. (2009). New sequence variants associated with bone mineral density. Nat. Genet. 41, 15–17.

Google Scholar

Tella, S. H., and Gallagher, J. C. (2014). Prevention and treatment of postmenopausal osteoporosis. J. Steroid Biochem. Mol. Biol. 142, 155–170.

Google Scholar

Tu, P., Duan, P., Zhang, R. S., Xu, D. B., Wang, Y., Wu, H. P., et al. (2015). Polymorphisms in genes in the RANKL/RANK/OPG pathway are associated with bone mineral density at different skeletal sites in post-menopausal women. Osteoporos. Int. 26, 179–185. doi: 10.1007/s00198-014-2854-7

PubMed Abstract | CrossRef Full Text | Google Scholar

Yang, C., Wang, Y., and Xu, H. (2017). Fluoride regulate osteoblastic transforming growth factor-beta1 signaling by mediating recycling of the type i receptor ALK5. PLoS One 12:e0170674. doi: 10.1371/journal.pone.0170674

PubMed Abstract | CrossRef Full Text | Google Scholar

Yang, T. L., Zhao, L. J., Liu, Y. J., Liu, J. F., Recker, R. R., and Deng, H. W. (2006). Genetic and environmental correlations of bone mineral density at different skeletal sites in females and males. Calcif. Tissue Int. 78, 212–217. doi: 10.1007/s00223-005-0267-5

PubMed Abstract | CrossRef Full Text | Google Scholar

Zheng, H. F., Forgetta, V., Hsu, Y. H., Estrada, K., Rosello-Diez, A., Leo, P. J., et al. (2015). Whole-genome sequencing identifies EN1 as a determinant of bone density and fracture. Nature 526, 112–117.

Google Scholar

Zhong, Z. A., Zahatnansky, J., Snider, J., Van Wieren, E., Diegel, C. R., and Williams, B. O. (2015). Wntless spatially regulates bone development through beta-catenin-dependent and independent mechanisms. Dev. Dyn. 244, 1347–1355. doi: 10.1002/dvdy.24316

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhu, D. L., Chen, X. F., Hu, W. X., Dong, S. S., Lu, B. J., Rong, Y., et al. (2018). Multiple functional variants at 13q14 risk locus for osteoporosis regulate RANKL expression through long-range super-enhancer. J. Bone Miner. Res. 33, 1335–1346. doi: 10.1002/jbmr.3419

PubMed Abstract | CrossRef Full Text | Google Scholar