Endochondral ossification pathway genes and postmenopausal osteoporosis: Association and specific allele related serum bone sialoprotein levels in Han Chinese

Osteoporosis is a systemic skeletal disorder characterized by reduced bone mineral density (BMD) and disrupted bone architecture, predisposing the patient to increased fracture risk. Evidence from early genetic epidemiological studies has indicated a major role for genetics in the development of osteoporosis and the variation in BMD. In this study, we focused on two key genes in the endochondral ossification pathway, IBSP and PTHLH. Over 9,000 postmenopausal Han Chinese women were recruited, and 54 SNPs were genotyped. Two significant SNPs within IBSP, rs1054627 and rs17013181, were associated with BMD and postmenopausal osteoporosis by the two-stage strategy, and rs17013181 was also significantly associated with serum IBSP levels. Moreover, one haplotype (rs12425376-rs10843047-rs42294) covering the 5’ end of PTHLH was associated with postmenopausal osteoporosis. Our results provide evidence for the association of these two key endochondral ossification pathway genes with BMD and osteoporosis in postmenopausal Han Chinese women. Combined with previous findings, we provide evidence that a particular SNP in IBSP has an allele-specific effect on mRNA levels, which would, in turn, reflect serum IBSP levels.

Scientific RepoRts | 5:16783 | DOI: 10.1038/srep16783 to be associated with BMD or OP 13 . Despite these attempts of association mapping based on common SNPs, the percentage of phenotypic variance explained by all these loci is still modest 13,14 . In addition, molecular and biological studies have lagged behind the genetic epidemiological studies. The biological mechanisms of over 30 BMD GWAS loci are still unclear 13 . Identifying specific susceptibility genes for BMD and OP can help reveal the underlying pathological mechanisms of OP, and in turn, inform the diagnosis, prevention, treatment strategy and novel drug development of this complex disorder.
Most bones form through the endochondral ossification, whereby the embryonic cartilaginous model gradually mineralizes and is replaced by bone 15 . As an essential process of skeletal development, it is systematically controlled by multiple regulatory factors. Genetic analysis incorporating the knowledge of biological pathways has identified gene clusters in the endochondral ossification pathway 16 . Several key factors of this pathway are associated with BMD or OP, including IBSP [17][18][19] , PTH 20 , RUNX2 21,22 , SOX6 23 and SP7 24 . Previous studies have shown that IBSP is significantly associated with BMD in Americans 17 and Australians 24 . Furthermore, unidirectional allele-specific regulation (ASE) of IBSP has been reported in human bone samples, with lower expression of the G allele compared to the A allele for SNP rs17013181 17 . However, it is still unclear whether the association of IBSP with BMD exists in a genetically independent Han Chinese population, and whether a similar allele-specific pattern identified at the RNA level would reflect serum protein levels.
To comprehensively investigate the susceptibility of genes involved in the endochondral ossification pathway, we examined the relationship between the two key genes IBSP and PTHLH and the diagnosis status of OP in over 9,000 postmenopausal women from the Han Chinese population. To determine whether the significant variants of OP are related to serum integrin-binding sialoprotein (IBSP, encoded by gene IBSP) levels, we measured IBSP serum levels in all of the samples.

Materials and Methods
Subjects. Two separate datasets were included in this study, and a two-stage approach was utilized for the discovery of single marker analyses. Subjects consisting of 1,425 women (aged 48-70 years) with primary postmenopausal osteoporosis and 3,557 healthy age-matched women (aged 48-70 years) were considered part of the discovery set and recruited from the Second Affiliated Hospital of Xi'an Jiaotong University and Xi'an Honghui Hospital. Subjects consisting of 1,039 women (aged 48-72 years) with primary postmenopausal osteoporosis and 3,032 healthy age-matched women (aged 48-72 years) were considered the replication set and enrolled from Taizhou Zhang's Orthopaedic Hospital. Two-stage subjects recruited in this study were random genetically unrelated Han Chinese postmenopausal women from the city of Xi'an in Shaanxi Province and Taizhou in Zhejiang province with no migration history within the previous three generations. None of the subjects had a history of taking medication to treat osteoporosis, and subjects with diseases or medications known to affect bone metabolism were excluded from the controls. Given the possibility of osteoporosis resulting from obesity and the overall level of Asian body mass index (BMI), subjects with a BMI ≥ 27 were excluded from the study. The demographic information of the study subjects is shown in Table 1. This study was performed in accordance with the   17,24 . Based on the above criteria, 54 SNPs were included in further analyses (Supplementary Table S1). Peripheral blood was drawn from a vein into a sterile tube containing ethylenediamine tetraacetic acid (EDTA). Genomic DNA was extracted from peripheral blood leukocytes according to the manufacturer's protocol (Genomic DNA kit, Axygen Scientific Inc., California, USA). DNA was stored at − 20 °C for SNP analysis. Genotyping was performed for all SNPs using the MassARRAY platform (Sequenom, San Diego, California, USA). Briefly, SNPs were genotyped using high-throughput, matrix-assisted laser desorption ionization-time-of-flight (MALDI-TOF) mass spectrometry. The resulting spectra were processed using Typer Analyzer software (Sequenom), and genotype data were generated from the samples. As the final genotype call rate of each SNP was greater than 98% and the overall genotyping call rate was 99.6%, the reliability of further statistical analysis was ensured.
Bone mineral density and serum bone turnover marker measurements. Dual-energy X-Ray absorptiometry (Lunar Expert 1313, Lunar Corp., USA) was used to measure BMD at the lumbar spine (L2-4) and femoral neck. Bone mineral density was determined according to standard Lunar protocols. Osteoporosis was defined according to the conventional World Health Organization (WHO) definition. Participants with a T score < − 1.0 SD were classified as having low bone mass, and those with a T score > − 1.0 SD were classified as normal. The anthropometric baseline data of all subjects were obtained through measurements and a questionnaire.
Peripheral venous blood was collected from each participant in the morning after a 12-hour fast according to the Declaration of Helsinki guidelines. Serum bone alkaline phosphatase (BALP) levels, osteocalcin (OST) levels, carboxyl-terminal telopeptide of type I collagen (CTX) levels, osteoprotegerin (OPG) levels, soluble receptor activator of nuclear factor-JB ligand (sRANKL) levels and bone sialoprotein levels (BSP) were measured using enzyme-linked immunosorbent assay (ELISA) kits (Westtang Biotech Inc., Shanghai, China). The minimal detection limit and the interassay and intra-assay coefficients of variation are presented for every marker in Supplemental Table S2. The summary of BMD and serum bone turnover marker measurements are shown in Table 1.

Statistical analysis. A Hardy Weinberg Equilibrium (HWE) test and single marker association anal-
ysis was performed for sample sets of both stages by Plink 26 . OP diagnosis was considered the main phenotype in single marker based analyses. BMD measurements of the lumbar spine and femoral neck were analyzed as two quantitative traits by fitting a linear model. Age and BMI were included as covariates when conducting single marker based analyses. Haplotype-based analyses were also performed using Plink 26 . Linkage disequilibrium (LD) blocks were constructed and plotted using Haploview 25 . Bonferroni correction was applied for both the discovery and replication stages. SNPs that passed the Bonferroni corrected P value threshold (0.05/54 ≈ 0.0009) in the discovery stage were genotyped in the replication sample sets. SNPs with significant hits by LD (r 2 > 0.1) in the discovery stage were included in the replication stage. In addition to the single marker based analyses, a conditioning model was fitted to test whether the significant SNPs identified in the single marker model were independent in the replication sample set. Furthermore, the association of serum IBSP level and genetic markers was analyzed by fitting a linear model using Plink 26 . All SNPs included in the replication stage were analyzed in combined sample sets. When fitting the linear model, age, BMI and sample site were considered covariates. A conditioning model was also fitted with the combined sample sets to test whether the effects of significant SNPs were independent. Bar plots were generated with R to visualize the association of SNP genotypes and quantitative traits 27 . Descriptive analyses were also conducted in R. Both Polyphen2 27 and SIFT 28 were used to predict the functional significance of significant non-synonymous SNPs identified in this study.

Results
Genetic association mapping of postmenopausal osteoporosis. The results of HWE and MAF are shown in Supplemental Table S1. All 54 SNPs passed the HWE test in the discovery stage. Two SNPs within the IBSP gene (rs1054627, P = 1.00 × 10 −4 ; rs17013181, P = 4.99 × 10 −5 ) passed the Bonferroni correction P value threshold in the discovery stage. We genotyped these 2 SNPs plus 11 additional SNPs in the same LD clusters in an independent analysis of the replication stage. The significance of both SNPs was successfully replicated in the second stage (rs1054627, P = 0.0001; rs17013181, P = 0.0003 with Bonferroni correction threshold of 0.05/13 ≈ 0.0038). After conditioning on SNP rs1054627, the significance of SNP rs17013181 was maintained (P = 0.0015). We did not identify any other SNPs with independent effects through single marker based analyses. The results of single marker based analyses of the 13 markers included in the replication stage are summarized in Table 2. We also performed association mapping analysis based on two BMD measurements (lumbar spine and femoral neck). The full results of association mapping for OP and BMD based on the discovery sample set are shown in Supplemental Table S3. The results of association mapping for BMD based on the replication sample set are presented in Supplemental Table S4. We confirmed the significance of both SNPs (rs1054627 and rs17013181) that we identified to be associated with OP in both stages in association analyses based on the BMD of the combined datasets (Fig. 1). As shown in Supplemental Figure S1 and S2, we constructed 18 LD blocks for the two genes (IBSP and PTHLH) based on the discovery sample set data. Haplotype analyses based on the discovery sample set data are summarized in Supplemental Table  S5. We identified some significant haplotypes in the IBSP gene (rs6828578-rs958848, P = 2.98 × 10 −6 ; rs4693878-rs1054627, P = 1.48 × 10 −7 ; rs13144371-rs17013181, P = 3.26 × 10 −7 ) and the PTHLH gene (rs12425376-rs10843047-rs42294, P = 5.75 × 10 −6 ).

Association of serum IBSP levels and IBSP gene SNPs.
We determined that several SNPs in the IBSP gene were significantly associated with serum IBSP levels in the combined sample set. However, after conditioning on the two most significant SNPs (rs17013181 and rs17013182), none of the other SNPs showed significance (  Table 2. Results of the single marker based association analysis for the 13 SNPs within IBSP genotyped in the replication stage. Significant SNPs were highlighted in bold. * Odds ratio and P values for analysis based on the sample of discovery stage. ** Odds ratio and P values for analysis based on the sample of replication stage. *** Odds ratio and P values for analysis based on the sample of replication stage while conditioning on SNP rs1054627. (r 2 = 0.83), we explored the relationship between serum IBSP levels and the combined genotype of the two SNPs (Fig. 2). After conditioning on the effects of rs17013181, we observed that the direction of the effects of rs17013182 was reversed (Table 3 and Fig. 2).  Table 3. Results of association analysis for serum IBSP level based on the data from combined sample set. Significant SNPs were highlighted in bold. * Linear regression coefficient and P values based on the data from combined sample set. ** Linear regression coefficient and P values after conditioning on SNP rs17013181. *** Linear regression coefficient and P values after conditioning on SNP rs17013181 and rs17013182.

Discussion
To our knowledge, this is the first large-scale genetic association study of IBSP and PTHLH genes with BMD, OP and serum IBSP levels in a Han Chinese population. Our results provide evidence for the association of two key endochondral ossification pathway genes with OP and BMD in postmenopausal Han Chinese women. Two SNPs with independent effects and some haplotypes in both IBSP and PTHLH genes were significantly associated with an OP diagnosis. We confirmed the effects of the two SNPs in analyses based on quantitative traits (BMD). Moreover, we also found evidence of an association between IBSP SNPs and serum IBSP levels.
Integrin-binding sialoprotein is an approximately 70 kDa acidic glycoprotein encoded by the human IBSP gene 29 . It belongs to the Small Integrin Binding Ligand N-linked Glycoprotein (SIBLING) gene family 30 . Previous studies have reported that IBSP is highly expressed in hypertrophic chondrocytes, which are at the center of endochondral ossification, and IBSP knockout mice have problems in bone resorption and mineralization 31 . In our study, we identified two non-synonymous SNPs within the IBSP gene (rs1701318 and rs1054627) that were significantly associated with the diagnosis of OP and variation in BMD. Conditioning analyses indicated that the effects of the two SNPs were independent. Previous reports have shown a significant association between SNP rs1054627 and BMD in populations of different ethnic groups. GWAS of Icelanders and replication samples of Europeans found this SNP to be significantly associated with femoral neck BMD 32 . It was also found to be significantly associated with femoral neck BMD in premenopausal European and African-American women 17 . In an extended genome-wide association study of bone mineral density based on 15,000 European individuals, rs1054627 in IBSP was significantly associated with bone density 18,19 . Among these previous studies, the directions of the effects were similar regardless of the ethnicity of the study samples. Although this SNP was predicted with limited functional significance by Polyphen2 (0.206) and SIFT (0.67), it might still alter protein structure and/or function because it changes the non-hydrophilic Gly residue to Glu, which is a hydrophilic residue 17 . Interestingly, this SNP was not associated with serum IBSP levels (P = 0.3134), which suggests that it did not affect IBSP expression level. Therefore, the effects of this SNP on OP diagnosis or BMD variation may solely come from the changes in protein structure. The other significant SNP, rs17013181, was also a non-synonymous SNP located in exon 7 of the IBSP gene. To our knowledge, this SNP have never been reported to be associated with OP or BMD diagnosis. In an earlier study 17 , researchers examined 52 human bone samples from the femoral neck during surgical hip replacement, and observed the allele-specific regulation of gene expression (ASE) for rs17013181. They reported that the minor allele G was expressed at lower levels than the ancestral allele A. In our study, we detected the association of rs17013181 with serum IBSP levels and observed a similar pattern ( Table 3). The minor allele G was associated with lower serum IBSP levels (β = − 0.6079, P = 4.51 × 10 −59 ). Our findings suggest that the ASE pattern of IBSP is reflected at the protein level. Similar to SNP rs1054627, this SNP was predicted to have limited functional significance by Polyphen2 (0) and SIFT (0.36). However, this SNP can still affect OP onset and development and BMD variation by altering IBSP expression levels. The rs17013182 SNP in strong LD with rs17013181 was also significantly associated with serum IBSP levels ( Table 3). Unlike rs17013181, this SNP only associated with serum IBSP levels, and not with the disease phenotype (Neither OP diagnosis nor BMD).
PTHLH is the protein-coding gene of parathyroid hormone-like related protein, which is a member of the parathyroid hormone family 33 . Parathyroid hormone (PTH) is a key regulator of calcium metabolism and contributes to skeletal development by regulating chondrocyte proliferation and differentiation during early bone growth 34 . Previous genetic studies have linked other PTH genes with age-related degenerative changes in lumbar spine 35 and BMD 20 . In this study, we identified a three-SNP haplotype (rs12425376-rs10843047-rs42294, P = 5.75 × 10 −6 ) that was significantly associated with the diagnosis of OP. The three SNPs were located near the 5′ end of the PTHLH gene. The functional relationship of the genomic region covered by this haplotype is unknown. However, this haplotype covers a transcription factor Chromatin Immunoprecipitation-Sequencing (ChIP-Seq) region. Several genes bind to this region, including EZH2, YY1, TBP, ESR1, MAX, REST, USF1, SIN3AK20, ZNF143, SP2, CEBPB and SIN3A 36 . Among these genes, ESR1 and EZH2 have been implicated in bone maturation 11,37,38 . However, further study is required to explain the significance of these association signals.
Our study has several strengths. The large sample size of our study (> 9,000) ensured the statistical power to detect modest effects of SNPs. In addition, we utilized a two-stage study design (discovery stage and replication stage). With this study design, SNPs that were identified as significant in both stages were unlikely to be false positives. Unlike previous studies that focused on the relationship between genetic polymorphisms and disease phenotypes, we measured serum IBSP levels to further investigate the potential effects of susceptibility SNPs on gene expression. However, several limitations should also be kept in mind. Possible population admixture and stratification could complicate our findings, although we have applied several strategies to reduce potential population stratification during the sample collection stage. Furthermore, our study consisted only of postmenopausal women, and from the perspective of epidemiology, this subject selection strategy impaired the representativeness of our study. We also failed to measure hormone levels in our subjects, which could confound our results. In addition, like most common SNP-based studies, we only included genetic variations with MAF ≥ 0.01 for association and therefore cannot detect the potential contribution of rare variants. Recent genetic epidemiology studies have indicated that rare variants play important roles in the onset and development of complex disorders 39 In addition to the genetic variations in its encoding gene, serum IBSP levels can be affected by several other biological factors that cannot be effectively examined or controlled in association analyses between serum IBSP level and SNPs in the IBSP gene. This could potentially weaken the convincingness of our findings.

Conclusions
In summary, our results provide further supportive evidence for the association of the IBSP and PTHLH genes with the diagnosis of OP and of the IBSP gene with variation of BMD in postmenopausal Han Chinese women. Combined with previous findings, we provided evidence that a particular SNP in IBSP has an allele-specific effect on the mRNA level, which would in turn, affect serum IBSP levels. Two SNPs in IBSP, rs1054627 and rs17013181, may confer risk of postmenopausal osteoporosis in Han Chinese women, and be useful in the informative assessment of the genetic risk for reduced BMD together with serum IBSP level. However, in view of multiple variants with small effects and the molecular basis of associations in relation to the complex network underlying bone remodeling and bone mass, further studies are required, including high density mapping and deep sequencing in different populations, to identify possible causal variants and determine the detailed mechanism of how these variants affect IBSP expression.