Common Variants in OPG Confer Risk to Bone Mineral Density Variation and Osteoporosis Fractures

Although many common variants have been identified for bone mineral density (BMD) and osteoporosis fractures, all the identified risk variants could only explain a small portion of heritability of BMD and osteoporosis fractures. OPG belongs to the tumor necrosis factor receptor superfamily, which plays a crucial role in bone remodeling and is thus a promising candidate gene of osteoporosis. Several studies have explored the association of OPG variants with BMD or osteoporosis fractures, however, the results remain inconsistent among different populations. In the study, we first assessed the relationship between OPG variants and BMD or osteoporosis fractures in our sample size (227 subjects with postmenopausal osteoporosis and 189 controls), and then performed a systematic meta-analysis. Among the nine SNPs genotyped, rs6469804 and rs2073618 showed significant associations with both BMD and osteoporotic fractures, while rs3102735 was only associated with BMD in our samples (P < 0.05). For meta-analyses, data for a total of 12 SNPs were pooled (4725 patients and 37804 controls), and five SNPs, including rs6993813, rs6469804, rs3134070, rs2073618 and rs3102734, showed association with osteoporosis fractures (P < 0.05). On light of the above analysis, we believe that OPG is one promising susceptibility gene of BMD or osteoporotic fractures.

Measurement of bone mineral density. BMD (grams per square centimeter, g/cm 2 ) at the lumbar spine (L1-L4) and total hip were measured by the dual-energy X-ray absorptiometry (DXA, GE Lunar Prodigy) for all subjects.
SNP selection and genotyping. SNPs across the OPG gene locus were selected based on the following three criteria. First, we chose those SNPs that have been reported to be assocaited with BMD or osteoporosis fractures in other populations. Second, part of tagging SNPs were retained. In order to choose tagging SNPs, the whole SNPs within OPG genomic region were downloaded from the Hapmap project (Chinese Han Populations) and the Haploview program (version 4.1) was applied to determine the tagging SNPs using the r 2 confidence interval (CI) algorithm 18 . The third was potential functional SNPs that might affect protein structure, mRNA expression, and so on. Finally, a total of 9 SNPs were selected for genotyping. The SNP Information was shown in Supplementary Table 1.
Genomic DNA was extracted from a total of 3 ml peripheral blood leukocytes using the standard phenol-chloroform method. Before genotyping SNPs, DNA was diluted to a final concentration of 30 ng/μl, and then all SNPs were genotyped by the TaqMan protocol as described in previous studies 19 . Genotyping assays were conducted on a CFX96 Real-Time PCR Detection System (Bio-Rad, Hercules, CA). To ensure accuracy, internal controls with known genotypes and negative controls with water were used. A total of 10% samples were randomly chosen for duplication to assess the genotyping error rate. The call rate for these SNPs was 97.4% on average and the genotype concordance was 100%.
Study searching and data collection for meta-analyses. Assocaition between OPG variants and fractures have been extensively studied before, and meta-analyses were performed acccording to previously reported procedures [20][21][22] . In brief, we searched eligible studies for this meta-analysis from PubMed (http://www. ncbi.nlm.nih.gov), SCOPUS (http://www.scopus.com), EMBASE (http://www.elsevier.com/online-tools/embase) and ISIWeb of Knowledge (http://apps.webofknowledge.com/) with the following searching items "(OPG or TNFRSF11B or TR1 or OCIF or PDB5) and fractures" in text. The last search was conducted in August 12 nd , 2016. Studies published in English were considered. We also screened the references of eligible studies for any other eligible studies.
Any studies with available genotyping data of OPG SNPs in patients with osteoporosis fractures (any location) as well as healthy controls were included. For each eligible study, two independent investigators (Sheng and Cai) extracted the following data: 1) first author and publication year; 2) fracture locations; 3) sample size and  softwares. Quantitative data was represented as mean ± SD. The chi-square (χ 2 ) goodness-of-fit test was applied to calculate the Hardy-Weinberg equilibrium (HWE) for each included SNP. Multiple regression and logistic regression analyses were carried out to evaluate the relationships between variables and bone mineral density as well as fractures risk, respectively. For meta-analyses, publication bias analysis, meta-analysis, and sensitivity analysis were investigated using the Stata 12.0 software. Either Egger regression test with a funnel plot or the Begg-Mazumdar test was applied to test potential publication bias. Cochran's χ 2 -based Q-statistic was calculated to assess the between-study heterogeneity. The pooled OR was determined by the Z-test using the random-or the fixed-effect model based on heterogeneity detection. For all test, P-value less than 0.05 was considered significant. Bonferroni correction was applied by multiplying the P-values by the number of tests whenever multiple tests were presented to surmount the inflation of P-value.

Results
OPG variants and BMD in our samples. A total of 9 SNPs were chosen for genotyping in our cohort, including 4 SNPs reported previously (rs6993813, rs6469804, rs3102735 and rs4355801), 3 tagging SNPs (rs7463176, rs1032128 and rs10955911) and 2 functional SNPs (rs2073617 and rs2073618). The detailed information, i.e., chromosome position, global and Chinese minor allele frequency (MAF), was listed in Supplementary  Table 1. The observed allelic frequency distributions in controls for all SNPs analyzed were in Hardy-Weinberg equilibrium (P > 0.05). Meanwhile, based on the following assumptions: P = 0.05, OR = 2.00 corresponding to a  Table 2. Association of 9 SNPs with BMD at Lumbar spine L1-L4 and total hip. a RAF, risk allele frequency. b P-value of significance was set at 0.006 by using the Bonferroni Correction. "moderate to high" effect, and a given MAF according to the 1000 Genomes project (www.1000genomes.org/), the present sample size revealed more than 80% statistical power for all SNPs (Supplementary Table 1). Table 2 showed the associations between the tested SNPs and BMD in particular locations. Three SNPs, namely rs6469804, rs3102735 and rs2073618, exhibited significant association with BMD in either Lumbar spine L1-L4 or total hip after adjusting for age, BMI, and age at menopause. Specifically, rs6469804 correlated with BMD with a P-value of 0.048 for Lumbar spine L1-L4 (β = −0.07, 95% CI = −0.23 to 0.00), and a P-value of 0.012 for total hip (β = −0.09, 95% CI = −0.25 to −0.04). While rs3102735 was only associated with total hip BMD (P = 0.006, β = −0.10, 95% CI = −0.28 to −0.07), and rs2073618 only with Lumbar spine L1-L4 (P = 0.015, β = −0.08, 95% CI = −0.30 to −0.04). After stringent Bonferroni correction, none of the significance remained (P-value for significance was set at 0.006). Table 3 showed the associations of each SNP with osteoporotic fractures. Corresponding to the association results of BMD, after adjusting for age, BMI, and age at menopause, rs6469804 (P = 0.026, OR = 0.656. 95% CI = 0.451-0.952) and rs2073618 (P = 0.039, OR = 0.721, 95% CI = 0.528-0.984) also showed significant associations with osteoporotic fractures, and the minor allele (rs6469804-G and rs2073618-G) showed protective effect for osteoporotic fracture onset. However, neither rs6469804 (P = 0.234) nor rs2073618 (P = 0.351) survived for Bonferroni correction (P-value for significance was set at 0.006).

Discussion
Like other genetically complex diseases, osteoporosis is the results of combinational effect of multiple genetic and environmental factors. Several studies have proved that genetic risk factors play a crucial role in the pathogenesis of osteoporosis with an estimated heritability of 55-85% 26,27 . By applying large-scale GWAS, many susceptibility loci conferring risk to BMD have been identified. For example, one recent published GWAS using European and east Asian cohorts identified 56 loci associated with BMD which surpassed genome-wide significance (P < 5E-08), and 14 SNPs among them showed significant association with fracture risk 7 . However, candidate association studies of BMD or osteoporotic fracture focusing on bone remodeling related genes have produced conflicting results in different populations 8,15 . Discordant results may be due to different populations or limited sample size used in previous studies. Success of genetic analysis of BMD and osteoporosis fracture will require either large sample sizes by meta-analysis.
In the present study, we first tested a total of 9 SNPs across the OPG locus which would influence the BMD and osteoporotic fracture in Chinese women, and found that 3 SNPs (rs6469804, rs3102735 and rs2073618) were associated with either spine or hip BMD after adjusting for other confounding facators, i.e., age, BMI, and age at menopause, further suggesting that common variants in OPG variation might contribute to BMD variation and even osteoporotic fracture. In order to validate our result, we then performed a comprehensive meta-analysis of the associations between OPG variation and osteoporotic fracture with the largest sample size to date (4725 patients and 37804 controls). In line with our expectations, a total of 5 SNPs showed significant association with osteoporotic fracture with the lowest P-value for rs3102734 (P = 6.3E-05).
With a key role in bone remodeling, OPG could bind to RANKL, a protein expressed on the osteoblast surface, thus blocking the interaction between RANKL and RANK on the osteoclast membrane. The OPG-RANKL interaction might exert antiresorptive effect in bone remodeling 9 . In vivo experiments also demonstrated that OPG knockout mice could develop severe osteoporosis 28 . OPG has been shown to be associated with osteoporosis fracture in various ethnicities. In a recent multi-ethnicity GWAS with more than 130000 cases and controls, authors found that SNPs in genes of the RANK-RANKL-OPG pathway despite showed the strongest BMD-associated signal 7 .
We could notice that not all variants that affect BMD are associated with osteoporosis fracture. As a matter of fact, only a small portion of BMD-related loci also show significant association with osteoporosis fracture. This might be not surprising when we take into consideration the following several reasons. First, although lower BMD is a prerequisite for fracture, only lower BMD combined with other risk factors could lead to osteoporotic fractures. Second, fracture is the extreme condition of osteoporosis.
Inconsistent with previous studies that enrolled men and women at the same time, we enrolled premenopausal patients only since it had been reported that estrogen deficiency might lead to bone loss and then confound the association studies 29 . This is one advantage of our study. There are, however, several limitations to the interpretation of our study. Firstly, although the combined sample size was large, most of the samples were from European districts, while Chinese samples account for only a small portion, which may contribute to type II errors of our results. Secondly, we focused only on the common variants in the OPG locus, and whether rare mutations in OPG also confer risk to BMD and even osteoporotic fractures remains largely unknown. Future study would focus more on rare variants of candidate genes through genome and candidate region sequencing. Thirdly, except for rs3134070 and rs3102734, the effect size for other risk variants were too small, and they might not represent causal variants, but only association signals which were in LD with other true causal variants. Fourthly, meta-analyses of the associations between OPG variants and BMD were not performed since the data necessary for meta-analyses was unavailable in most of studies. In this case, the meta-analytic results might be inaccurate since sample variance couldn't be deduced from the original manuscript.

Conclusion
In light of the above analysis, we believe that OPG is one promising susceptibility gene of BMD or osteoporotic fractures. However, whether rare variants also confer risk to BMD or osteoporotic fractures remain to be identified.