Genetic risk score based on the prevalence of vertebral fracture in Japanese women with osteoporosis.

A genetic risk score (GRS) was developed for predicting fracture risk based on the prevalence of vertebral fractures in 441 Japanese females with osteoporosis. A total of 979 (858 nonsynonymous and 121 silent) single-nucleotide polymorphisms (SNPs) located in 74 osteoporosis-susceptibility genes were genotyped and evaluated for their association with fracture prevalence. Four SNPs (protein kinase domain containing, cytoplasmic [PKDCC; rs4952590], CDK5-regulatory subunit-associated protein 1-like 1 [CDKAL1; rs4712556], wingless-type MMTV-integration site family member 16 [WNT16; rs2707466], and G-patch domain-containing gene 1 [GPATCH1; rs10416265]) showed a significant association (p < 0.05) with the fracture, in which the minor allele of the former two SNPs was the protective allele and that of the latter two SNPs was the risk allele. Applying a dominant-genetic model, we allotted - 1 point each to the protective-allele carriers and 1 point each to the risk-allele carriers, and GRS values were calculated as the sum of the points. The receiver-operating characteristic curves showed that GRS adequately predicted vertebral fracture. For the model predicted by the GRS with and without the effect of age, areas under the curves were 0.788 (95% confidence interval [CI]: 0.736-0.840) and 0.667 (95% CI: 0.599-0.735), respectively. Multiple logistic regression analysis revealed that the odds ratio for the association between fracture prevalence and GRS was 3.27 (95% CI: 1.36-7.87, p = 0.008) for scores of - 1 to 0 (n = 303) and 12.12 (95% CI: 4.19-35.07, p < 0.001) for scores of 1 to 2 (n = 35) relative to a score of - 2 (n = 103). The GRS based on the four SNPs could help identify at-risk individuals and enable implementation of preventive measures for vertebral fracture.

Introduction

Osteoporosis is among the most common skeletal diseases and affects > 200 million individuals worldwide, a figure that continues to increase as populations in developed countries live longer than previous generations (Reginster and Burlet, 2006). Osteoporosis is clinically characterized by reduced bone mass and compromised bone strength, which leads to an increased risk of fracture (Soen et al., 2013).

Fragility fractures, such as vertebral and femoral fractures, are among the most serious complications in elderly patients with osteoporosis. Several important factors, including age, past history of fragility fractures, family history of femoral fracture, bone mineral density (BMD), and history of falls, increase the risk of fracture in a clinical setting. In addition to these factors, genetic variations also determine predisposition to low-trauma fractures as demonstrated by genetic–epidemiological studies (Peacock et al., 2002, Ralston and Uitterlinden, 2010). Recent large-scale meta-analyses of genome-wide association studies (GWAS) identified a number of single-nucleotide polymorphisms (SNPs) associated with low BMD or increased fracture risk (Styrkarsdottir et al., 2008, Rivadeneira et al., 2009, Estrada et al., 2012).

Risk scores have been developed to predict the risk of coronary heart disease (ERICA Research Group, 1991, Tunstall-Pedoe, 1991), diabetes mellitus (Lindström and Tuomilehto, 2003), and dementia (Ngandu et al., 2006), and typically encompass multiple factors that affect disease onset or progression. A risk score can therefore improve the ability to predict common polygenic diseases, such as osteoporosis, by including multiple SNP profiles. Accordingly, we previously performed a study to develop a genetic risk sore (GRS) for predicting lifetime femoral fracture risk using data from consecutive, elderly Japanese autopsy cases at a community-based geriatric hospital (Zhou et al., 2015). The aim of this study was to develop a GRS for predicting vertebral fracture risk in Japanese women with osteoporosis using data registered in Biobank Japan (Nakamura, 2007).

Methods

Subjects

The osteoporosis-case subjects were collected under the support of the BioBank Japan Projects (Nakamura, 2007), and all participants provided written informed consent as approved by the ethics committees of the BioBank Japan Project (Nakamura, 2007) and the University of Tokyo. Osteoporosis was diagnosed based on the Japanese diagnostic criteria for primary osteoporosis (Soen et al., 2013). Patients with malignant neoplasms, liver cirrhosis, nephrotic syndrome, diabetes mellitus, rheumatoid arthritis, cerebral infarction, chronic obstructive pulmonary disease, hyperthyroidism, renal failure, and history of steroid-drug use were excluded from the assessment. Finally, 441 unrelated females with a mean age of 69.6 years were selected for this study. The prevalence of morphological vertebral fracture in all study subjects was determined by examination of lateral thoracolumbar (T4–L4) radiographs. The assessment of vertebral fracture was made in accordance with the semi-quantitative method (Genant et al., 1993), and a vertebral fracture was defined as a deformity of more than grade 1 in any of the measured vertebrae. Of the 441 subjects, 72 individuals sustained vertebral fractures, with mean age and age distributions (standard deviation and min–max, respectively) of 74.5 years (7.1 and 53–88) for subjects with fractures and 68.0 years (8.2 and 28–88) for those without.

SNP selection and genotyping

A large-scale meta-analysis of previous GWASs identified 56 BMD loci and revealed 14 loci associated with fracture risk (Estrada et al., 2012). To select SNPs for this study, those within or close to the 56 BMD loci were evaluated, as well as those on the Illumina HumanExome BeadChip (Grove et al., 2013) (Illumina, Inc., San Diego, CA, USA). A total of 979 (858 nonsynonymous and 121 silent) SNPs in 74 genes were identified (Supplementary Table 1) and evaluated for their association with the incidence of vertebral fracture among the 441 cases. The genotyping data for the 979 SNPs of the study subjects were provided from the Biobank Japan genotyping database generated using Illumina OmniExpressExome BeadChip version 1.2 (Illumina, Inc.) with call rates of > 0.99 during the process of genotyping.

Calculation of GRS

GRS was calculated as reported previously (Zhou et al., 2015). In this study, we applied a dominant-genetic model and allotted − 1 point each to the protective-allele carriers and 1 point each to the risk-allele carriers, and unweighted GRS values were calculated as the sum of the points. We also standardized scores using coefficients obtained from the logistic regression analyses (weighted GRS) to ensure that the lowest absolute value of the coefficient was assigned a value of 1 (Zhou et al., 2015). The association of the GRS with vertebral fracture was evaluated by multiple logistic regression analysis.

Statistical analysis

All statistical analyses were carried out using PLINK 1.07 software (http://pngu.mgh.harvard.edu/purcell/plink/) (Purcell et al., 2007) or SPSS for Windows version 20 (SPSS Inc., Chicago, IL, USA). Allelic frequencies of the selected SNPs were calculated using a gene-counting method. Hardy-Weinberg equilibrium for each SNP was assessed by the χ2 test. The Cochran-Armitage proportion trend test was used to identify changes in fracture incidence with respect to the number of risk or protective alleles. Multiple logistic regression analysis, including age, genotypes of each SNP, or GRS as independent variables, was performed to estimate the odds ratio (OR) and 95% confidence interval (CI) for the association with risk of vertebral fracture. Receiver-operating characteristic (ROC) curves were plotted, and the area under the curve (AUC) was calculated to assess the discriminative power of the GRS models (Ngandu et al., 2006). All reported p-values are two-sided, with p < 0.05 regarded as statistically significant.

Results

After genotyping 979 SNPs located in 74 previously reported osteoporosis-susceptibility genes (Estrada et al., 2012), SNPs that met the following criteria were selected: 1) a statistically significant association (p < 0.05) with the prevalence of vertebral fracture according to the Cochran-Armitage trend test, and (2) a minor-allele frequency of > 0.01 in the study population. Finally, four SNPs [protein kinase domain containing, cytoplasmic (PKDCC; rs4952590), CDK5-regulatory subunit-associated protein 1-like 1 (CDKAL1; rs4712556), wingless-type MMTV-integration site family member 16 (WNT16; rs2707466), and G-patch domain-containing gene 1 (GPATCH1; rs10416265)] were selected for this study (Table 1). Their allele and genotype frequencies were in Hardy-Weinberg equilibrium, and the prevalence of vertebral fracture was calculated for each genotype to identify risk or protective alleles. As shown in Table 1, the T allele of rs4952590 and the A allele of rs4712556 significantly protected vertebral fracture (protective allele), while the T allele of rs2707466 and the A allele of rs10416265 contributed significantly to fracture morbidity (risk allele).

Table 1

SNPs selected to calculate GRS for 441 female subjects with osteoporosis and the prevalence of vertebral fracture for each genotype.

Gene	SNP	Genotype	n	Vertebral fracture, n (%)	p
PKDCC	rs4952590	CC	144	32 (22.2)	0.042⁎
		CT	221	30 (13.6)
		TT	76	10 (13.2)
CDKAL1	rs4712556	GG	150	35 (23.2)	0.035⁎
		AG	219	26 (11.9)
		AA	72	11 (15.3)
WNT16	rs2707466	CC	339	47 (13.9)	0.012⁎
		CT	98	24 (24.5)
		TT	4	1 (25.0)
GPATCH1	rs10416265	GG	302	42 (13.9)	0.042⁎
		AG	127	27 (21.3)
		AA	12	3 (25.0)

CDKAL1, CDK5-regulatory subunit-associated protein 1-like 1; GPATCH1, G-patch domain-containing gene 1; GRS, genetic risk score; PKDCC, protein kinase domain containing, cytoplasmic; SNP, single-nucleotide polymorphism; WNT16, wingless-type MMTV integration-site family member 16.

Significant (p < 0.05) according to the Cochran-Armitage trend test.

The independent association of each minor allele with vertebral fracture was evaluated by multiple logistic regression analysis. A minor allele-dominant genetic model was used for these four SNPs, because the fracture prevalence was quite similar between heterozygous and homozygous carriers (Table 1). As shown in Table 2, T-allele carriers of rs4952590 and A-allele carriers of rs4712556 showed significantly lower ORs for risk of vertebral fracture, whereas T-allele carriers of rs2707466 and A-allele carriers of rs10416265 showed significantly higher ORs for the risk. In order to calculate GRS values, we allotted − 1 point each to the protective-allele carriers and 1 point each to the risk-allele carriers (unweighted-risk score). We also standardized scores using coefficients obtained from logistic regression analyses to ensure that the lowest absolute value of the coefficient was assigned a value of 1 (weighted-risk score) and found that the unweighted- and weighted-risk scores were quite similar (Table 2). Therefore, the unweighted-risk score was used for subsequent analyses.

Table 2

Multiple logistic regression analysis of the association between vertebral fracture prevalence and genotype.

	Coefficient	p	OR	95% CI	Unweighted-risk score	Weighted-risk score
rs4952590, CC	0 (reference)		1		0	0
rs4952590, CT/TT	− 2.024	0.043⁎	0.57	0.33–0.98	− 1	− 1
rs4712556, GG	0 (reference)		1		0	0
rs4712556, AG/AA	− 2.611	0.009⁎	0.49	0.28–0.84	− 1	− 1.29
rs2707466, CC	0 (reference)		1		0	0
rs2707466, CT/TT	2.092	0.036⁎	1.86	1.04–3.32	1	1.03
rs10416265, GG	0 (reference)		1		0	0
rs10416265, AG/AA	2.517	0.012⁎	2.05	1.17–3.59	1	1.24

CI, confidence interval; OR, odds ratio.

Significant (p < 0.05) according to multiple logistic regression analysis after adjusting for age.

Table 3 shows the prevalence of vertebral fracture based on the unweighted GRS, which was calculated as the sum of the unweighted-risk scores. The OR for the risk of vertebral fracture tended to increase with the number of the GRS, but was apparently influenced by the uneven distribution of the number of subjects. When the subjects were divided into three groups according to the GRS, namely those with GRS of − 2, − 1 to 0, and 1 to 2, the ORs were 3.27 and 12.12 for the groups with GRS of − 1 to 0 and 1 to 2, respectively, relative to the group with GRS of − 2 (Table 3). The ROC curves showed that the GRS adequately predicted vertebral fracture (Fig. 1). For the model predicted by the GRS with and without the effect of age, AUCs were 0.788 (95% CI: 0.736–0.840) and 0.667 (95% CI: 0.599–0.735), respectively.

Table 3

Risk of vertebral fracture based on unweighted GRS.

GRS	n	Vertebral fracturen (%)	OR (95% CI)	p
− 2	103	6 (5.8)	1 (reference)
− 1	180	26 (14.4)	2.73 (1.08–6.87)	0.033⁎
0	123	25 (20.3)	4.12 (1.62–10.50)	0.003⁎
1	31	14 (45.2)	13.31 (4.49–39.46)	< 0.001⁎
2	4	1 (25.0)	5.39 (0.48–59.92)	0.171
− 1 to 0	303	51 (16.8)	3.27 (1.36–7.87)	0.008⁎
1 to 2	35	15 (42.9)	12.12 (4.19–35.07)	< 0.001⁎

CI, confidence interval; GRS, genetic risk score; OR, odds ratio.

Significant (p < 0.05) according to multiple logistic regression analysis after adjusting for age.

Fig. 1

ROC curves of unweighted GRS with and without the effect of age, for predicting vertebral fracture risk.

Discussion

The main finding of this study was that the prevalence of vertebral fracture in 441 female patients with osteoporosis increased in proportion to the increase in the GRS, which was calculated based on the risk or protective-allele profiles of four SNPs (PKDCC rs4952590, CDKAL1 rs4712556, WNT16 rs2707466, and GPATCH1 rs10416265). The ROC curves showed that GRS adequately predicted the fracture risk, with AUCs with and without the effect of age of 0.788 (95% CI: 0.736–0.840) and 0.667 (95% CI: 0.599–0.735), respectively. The prevalence of vertebral fracture in the group with GRS of − 1 to 0 was 16.8% (Table 3), which was almost identical to that of all subjects enrolled in this study (72/441, 16.3%), indicating that the group with GRS of − 1 to 0 (n = 303) represented the ordinary risk group, whereas the groups with GRS of − 2 (n = 103) and 1 to 2 (n = 35) represented low- and high-risk groups, respectively.

Here, the T allele of WNT16 rs2707466 was found to be a risk allele for vertebral fracture. WNT proteins belong to a family of secreted cysteine-rich glycoproteins that transmit signals through both the WNT–β-catenin pathway, also termed the canonical WNT pathway, and noncanonical WNT pathways (Kohn and Moon, 2005). Within the components of WNT signaling, the gene encoding WNT16, one of the 19 WNT ligands in the human genome, is strongly associated with specific bone traits, such as cortical bone thickness, cortical porosity, and fracture risk (Baron and Kneissel, 2013, Movérare-Skrtic et al., 2014). Zheng et al. (2012) identified a novel missense SNP (C > T; Thr > Ile; rs2707466) located in the WNT16 gene and associated with cortical bone thickness by performing two separate GWAS meta-analyses in three cohorts comprising 5878 European subjects. Niu et al. (2016) conducted a three-stage meta-analysis targeting phosphorylation-related SNPs for femoral neck-BMD, total hip-BMD, and lumbar spine-BMD phenotypes, and found that WNT16 rs2707466 was associated with BMD phenotypes in each respective stage and in three stages combined, achieving genome-wide significance for both femoral neck- and total hip-BMD. In silico analyses predicted that rs2707466 directly abolishes a phosphorylation site, which could cause a deleterious effect on the WNT16 protein (Niu et al., 2016). WNT16 rs2707466 also influences heel-bone properties in a population of young adults as measured by quantitative-ultrasound techniques, which revealed aspects of bone fragility distinct from BMD (Correa-Rodríguez et al., 2016).

GPATCH1 is a gene of unknown function; however, the SNP selected for this study, GPATCH1 rs10416265, is a genetic determinant of heel-bone properties as determined by broadband-ultrasound attenuation and velocity of sound according to a GWAS meta-analysis (Moayyeri et al., 2014). PKDCC encodes a protein kinase belonging to a category of secretory pathway kinases that phosphorylate proteins and proteoglycans in the secretory pathway and appear to regulate various extracellular processes (Sreelatha et al., 2015); Imuta et al. (2009) used gene-knockout techniques to show that a protein-kinase gene, Pkdcc (AW548124), is required for longitudinal bone growth by promoting the appropriate differentiation of chondrocytes. A CDKAL1 gene variant, rs7756992, which is different from CDKAL1 rs4712556 selected in this study, influences insulin response and risk of type 2 diabetes (Steinthorsdottir et al., 2007). However, given the limited knowledge of the roles of these genes in bone biology, it is difficult to completely exclude the possibility that they are not causal and are in linkage disequilibrium with other bona fide osteoporosis-susceptibility genes.

Tran et al. (2011) first suggested that genetic profiling could enhance the predictive accuracy of fracture prognosis based on the clinical data obtained from the Dubbo Osteoporosis Epidemiology Study as well as simulated genetic data of 50 independent genes with allele frequencies ranging from 0.01 to 0.60 and relative risks ranging from 10.1 to 3.0. A recent large-scale meta-analysis identified 63 autosomal SNPs associated with BMD, of which 16 were also associated with fracture risk (Estrada et al., 2012). Estrada et al. (2012) evaluated the combined effect of the 63 BMD-associated SNPs to predict the risk for osteoporosis and fracture based on genotyping data obtained in the Prospective Epidemiological Risk Factor study, a prospective study of postmenopausal Danish women (Bagger et al., 2007). This study represented an independent-validation setting and was excluded from the overall meta-analysis for this reason. Despite serving as robust proof of the relationship between BMD-decreasing alleles and the risk of osteoporosis and fracture, prediction ability was modest. ROC analysis showed a significant, but relatively small discrimination ability of the genetic score alone, with AUCs of 0.59 and 0.57 for osteoporosis and fracture, respectively. Eriksson et al. (2015) also developed two GRSs, GRS63 and GRS16, based on the 63 BMD-associated and 16 fracture-associated SNPs, respectively, in order to determine the clinical usefulness of these GRSs for the prediction of BMD and fracture risk in elderly subjects. They studied two male and one female large prospective cohort of older subjects and found that GRS63 was associated with BMD and both GRS63 and GRS16 were associated with fractures. However, after BMD adjustment, the effect sizes for these associations were substantially reduced, and they concluded that, when BMD is known, the clinical utility of the two GRSs for fracture prediction is limited in elderly subjects. Lee et al. (2013) developed a GRS including 21 SNPs in 19 osteoporosis-susceptibility genes, and demonstrated that adding the GRS to the prediction model consisting of clinical risk factors and BMD could improve its predictive ability for non-vertebral fracture in 1229 unrelated Korean postmenopausal women. Lee et al. (2016) also calculated the Korean-specific GRS from 35 SNPs associated with osteoporosis-related traits (GRS35), and found that integration of the GRS35 into the current model further improved its predictability for future osteoporotic fracture occurrence in a 6-year follow-up observational study. The 979 SNPs genotyped in the present study were located in the 74 genes previously reported as osteoporosis-susceptibility genes (Estrada et al., 2012). However, it should be noted that most of the 979 SNPs were not identical to the originally reported marker SNPs used in the GWAS, which were located mainly in noncoding regions. In fact, the four SNPs selected in this study were not found in the 63 BMD-associated SNPs. Therefore, our GRS could not be directly compared with previous GRSs. Furthermore, given that it was a retrospective, patient-based study, the predictive value of our GRS for vertebral fracture should also be investigated in a prospective, population-based study.

Among the 979 SNPs tested, 133 SNPs were found to have a minor-allele frequency of > 0.01 in the present study population, and thus the p-values for the 4 SNPs were higher than the threshold significance level of Bonferroni correction for multiple testing (α = 0.05/133 SNPs = 0.000376). We also calculated the q-value (false discovery rate) for each SNP (Benjamini and Hochberg, 1995): 0.37 for rs10416265 (GPATCH1), 0.37 for rs2707466 (WNT16), 0.48 for rs4712556 (CDKAL1), and 0.48 for rs4952590 (PKDCC), which were not statistically significant. Therefore, we could not completely rule out the possibility of false positivity. However, it is also possible that the sample size of 441 subjects in the present study was relatively small and not sufficient to eliminate statistical ambiguity. Further studies including larger study samples are needed to elucidate the possibilities.

Vertebral fracture is among the most serious complications in elderly osteoporosis patients. A few important factors, including age, past history of fragility fractures, family history of femoral fracture, BMD, and history of falls, increase the risk for vertebral fracture in a clinical setting. Information regarding these factors along with a GRS based on risk-allele profiles of the four SNPs could therefore help to identify at-risk individuals to ensure that measures for preventing vertebral fracture can be implemented.

The following is the supplementary data related to this article.

Supplementary Table 1

Nine hundred and seventy nine SNPs located in the 74 osteoporosis susceptibility genes equipped on the Illumina HumanExome BeadChip.

Conflicts of interest

Heying Zhou, Seijiro Mori, Tatsuro Ishizaki, Atsushi Takahashi, Koichi Matsuda, Yukihiro Koretsune, Shiro Minami, Masahiko Higashiyama, Shinji Imai, Kozo Yoshimori, Minoru Doita, Akira Yamada, Satoshi Nagayama, Kazuo Kaneko, Satoshi Asai, Masaki Shiono, Michiaki Kubo, and Hideki Ito declare that they have no conflicts of interest.