Evaluation of WISP1 as a candidate gene for bone mineral density in the Old Order Amish

Wnt1-inducible signaling pathway protein-1 (WISP1) is a novel target of the Wnt pathway for modulating osteogenesis and improving bone strength. However, it is not clear if genetic variants in the WISP1 region are associated with bone mineral density (BMD) in human. The aim of this study is to investigate the role of genetic variation in WISP1 gene as a determinant of BMD in 1,510 Old Order Amish (OOA). We performed regional association analysis of 58 tag variants within 5 kb upstream and downstream to WISP1 with BMD and found 5 variants that were associated with BMD at multiple skeletal sites (P values from 2.89 × 10−6 to 1.62 × 10−2), with some significant associations even after adjustment for multiple comparisons. To replicate these results in an independent dataset, we performed a look-up of BMD associations with these variants in European ancestry subjects from the large GEFOS Consortium and observed the nominal associations of two of these variants with BMD (P values: 0.031 to 0.048). In conclusion, we have demonstrated that genetic variants surrounding WISP1 are associated with BMD at multiple skeletal sites in the OOA, thus influencing osteoporosis risk. These results support a role for the WISP1 gene on influencing variation in BMD.

Heritability of the BMD at multiple skeletal sites in the OOA. The heritability of each of the densitometric phenotypes was shown in Table 2 and was based on the most parsimonious model of variance component analysis for each phenotype, including only significant sources of variation. All of the results were statistically significant (p < 0.05). In the whole group consisting of males and females, the heritabilities of the BMD measurements in specific sites were high with variations from 0.66 (lumbar spine BMD) to 0.58 (hip femoral neck BMD). In the sex-stratified model, we found that heritabilities of BMD at multiple skeleton sites were generally greater in females than in males. The heritability of hip intertrochanter BMD (h 2 = 0.56) was the highest in the male group, whereas the heritability of lumbar spine BMD (h 2 = 0.67) was the highest in female group.
Replication study in GEFOS. We further looked for associations of the 5 lead SNPs in the published meta-analysis data in Caucasian 20 . We found that one SNP (rs11778573) were nominally associated with BMD at both femoral neck (p = 4.82 × 10 −2 ) and lumbar spine (p = 3.20 × 10 −2 ), and another SNP (rs116873248) showed suggestive level of association with BMD at lumbar spine (p = 3.10 × 10 −2 ) ( Table 3). We compared the minor allele frequencies of these 5 polymorphisms in the OOA with minor allele frequencies in the GEFOS Caucasian population and found that OOA minor allele frequencies varied little from GEFOS Caucasian population allele frequencies (Table 3).  Conditional analyses. To further determine independent variants, we performed a conditional analysis using the top variant rs72731533 as a covariate and no signal in this region remained (p > 0.01) (Fig. 2). This result indicated that only one association signal in this region was associated with BMD at multiple skeletal sites. We performed linkage disequilibrium (LD) analysis for the five variants associated with BMD using Haploview and found that rs35513885 (exon 4, A118S) was in high LD with rs72731533 (r 2 = 0.70, Fig. 3). Using rs35513885 as a covariate, conditional association analysis showed that rs72731533 was still associated with BMD traits (  Fig. 4D, rs72731533 is situated near the enhancer elements (H3K4Me1 and H3K27Ac marks) and also near the promoter elements (H3K4Me3 and H3K9Ac marks) on the BMDMSC (Bone Marrow Derived Mesenchymal Stem Cell Cultured Cells) and cfBMDMSC (Chondrocytes from Bone Marrow Derived Mesenchymal Stem Cell Cultured Cells) from ENCODE, which indicated that rs72731533 was probably involved in the regulation of gene expression. The SNP, rs72731533, also fell into a DNase Hypersensitivity site in primary osteoblasts (Fig. 4E,F). The accessible chromatin zone is functionally related to transcriptional activity, since this remodeled state is necessary for the binding of proteins such as transcription factors. We further examined whether the 5 significant association SNPs influence gene expression using public databases such as the Genotype-Tissue Expression (GTEx) project and ExSNP. However, we did not find significant eQTL for these SNPs that may be due to lack of bone tissue information in the databases. We further conduct  cis-eQTL analyses on the SNPs in WISP1 gene region with transcripts in primary osteoblasts (obtained from bone biopsies) 21 . We found that a SNP, rs144161059, in high LD (linkage disequilibrium) with rs72731533 (D' = 1) was significantly associated with WISP1 gene expression (P = 2.17 × 10 −3 ). The SNP, rs144161059, with low minor allele frequency (MAF = 0.003 in 1000 Genomes) failed to be imputed in this study.

Discussion
The Wnt signaling pathway is one of the most important signaling pathways in bone regulation because of its essential role in development, particularly in tissue specification and in cellular migration 22 . This signaling pathway influences all types of bone cells (osteoblasts, osteoclasts and osteocytes) and has showed to be important in skeletal development, maintenance of skeletal homeostasis and in bone remodeling 23 . Recently, several independent studies, with the goal to detect candidate genes underlying osteoporosis, revealed that many genes in the Wnt signal pathway are associated with lumbar spine, hip femoral neck and whole body BMD, bone strength, cortical bone thickness, and fracture risk 10,24 . Within the components of the Wnt signaling pathway, the gene coding for WISP1 has been found as a novel target for modulating osteogenesis and improving bone strength. The importance of WISP1 gene in the regulation of BMD and bone strength had been also confirmed by WISP1 knockout (WISP1 −/− ) mice 14 and human WISP1 transgenic mice 15 . In this study, we have investigated whether the variants in the WISP1 gene are associated with BMD at multiple skeletal sites in the OOA. Our finding is the first report that shows the significant association of the variant surrounding WISP1 gene on BMD in the OOA. The heritability of BMD is significantly high (h 2 >0.5) reported by several studies [6,17,27,30]. We calculated heritability of BMD at multiple skeletal sites in 1,510 OOA subjects and compared the difference of heritability between males and females. The heritability of lumbar spine BMD (h 2 = 0.66) was the highest and the heritability of hip femoral neck BMD (h 2 = 0.58) was the lowest in the whole group of males and females combined suggesting a greater genetic determinant of BMD in the lumbar spine than in proximal femur that was similar to previous studies 25 . Variance in the BMD heritability of different skeletal sites is possibly due to dissimilar external forces placed on certain bones of the skeleton. Previous studies had successfully noted site-specific variations of BMD heritability that are gender-dependent 26 . In this study, we observed that both female and male groups have strong genetic correlations of BMD and we found that heritability of BMD was partly different between females (h 2 :0.59-0.69) and males (h 2 :0.53-0.56). It was similar to our previous result that the heritability in BMD was larger in women than in men 27,28 . However, our results do not support some previous studies that males tend to have higher heritability values than females 29,30 . Amish males perform farm work and Amish females perform housekeeping work. There is higher physical activity in Amish males than in Amish females that may contribute the difference of heritability in BMD. Tse and colleagues observed that the high degree of BMD heritability may reflect a preservation of genetic influence in the relative absence of external forces 26 . Our results revealed that 53-69% of variance in BMD of OOA is contributed by genetic factors.
WISP1 is a member of the CCN family growth factors and its variants have been confirmed to be associated with various pathologies. Tao J. et al. 17 found that the WISP1 rs16893344 variant allele (T) was associated with a significantly increased risk of myocardial infarction. Chen J. et al. 19 found that WISP1 rs16893344, rs2977530, rs2977537, and rs62514004 (P < 0.05) polymorphisms were related to susceptibility of lung cancer; and WISP1 rs11778573, rs16893344, rs2977536, rs2977549 and rs62514004 polymorphisms were significantly associated with platinum-based chemotherapy response in lung cancer patients. Interestingly, Tomohiko et al. 16 evaluated the association of a SNP (rs2929970) in the WISP1 3′-UTR region with the presence of osteophytes, endplate sclerosis, and narrowing of disc spaces for spinal osteoarthritis in 304 postmenopausal Japanese women and found strong associations of rs2929970 with endplate sclerosis. Several GWAS identified BMD SNPs are nominally associated with prevalent radiographic knee osteoarthritis (OA) 31 . The previous studies suggested that the loci, associated with osteoarthritis, might be also association with BMD. In the present study, we examined the association between polymorphisms around the WISP1 gene and BMD at the total hip, femoral neck, hip intertrochanter, hip trochanter and lumbar spine to investigate a possible contribution of WISP1 to human bone metabolism.
In this study, we identified 58 tag genetic variations with a minor allele frequency (MAF) of at least 1% surrounding the WISP1 gene (+/−5 kb) and found 5 variants significantly associated with all BMD traits. The top associated variant was a SNP (rs72731533) that was located in intron 2 of the WISP1 gene. We further investigated whether the 5 SNPs are independently associated with BMD and found that only one association signal in this region is associated with BMD at multiple skeletal sites. A previous study reported that an intergenetic SNP, rs7839059, located on chromosome 8q24.12 near to WISP1 gene was significantly associated with cortical vBMD (P = 1.2 × 10 −10 ) 32 . We found that rs7839059 was associated with all the phenotypes (total hip BMD, p = 5.32 × 10 −5 ; hip femoral neck BMD, p = 2.88 × 10 −4 ; hip intertrochanter BMD, p = 1.52 × 10 −4 ; hip trochanter BMD, p = 6.06 × 10 −5 and lumbar spine BMD, p = 3.53 × 10 −2 ). We performed a conditional analysis that showed rs7839059 was independent of the associated signal in WIPS1 gene (Table 5). To replicate our results in larger sample sets, we checked the published meta-analysis data in Caucasians 20 and found that two SNPs (rs11778573 and rs116873248) were nominally associated with BMD at multiple skeleton sites (p = 0.031-0.0482). We think the following two reasons may cause nominal replications. First, we used chips for genotyping that just included common variants. The significantly associated variants may be partially link to causal variants. Second, the  significantly associated variants may be population specific. The most significant SNP, rs72731533, was involved in regulatory elements (such as enhancer and promoter around WISP1 gene) in both MSCs and OPCs. We did not find significant eQTL for rs72731533. However, the SNP, rs144161059, in high LD with rs72731533 was significantly associated with WISP1 gene expression in primary osteoblasts, although, we did not imputed rs144161059 in this study. Thus, a denser fine-mapping study using sequencing of the WISP1 locus will provide a better resolution to identify potential causal variant(s) and will be helpful for future functional validation. Those reports suggested that variants around the WISP1 region were actively involved in regulation of multiple phenotypes. This combined evidence suggests that polymorphisms around the WISP1 are associated with BMD at multiple skeletal sites.
In conclusion, we performed a regional analysis for 5 kb upstream and downstream WISP1 with specific BMD adjusted for study, sex and age in 1,510 OOA. We confirmed that genetic variation at the WISP1 locus is significantly associated with BMD at multiple skeletal sites. These results identify that genetic variants in WISP1 gene region are associated with BMD levels. Bioinformatics analyses suggest that this feasible association is partly caused by regulatory effects on the enhancer or promoter of WISP1 gene. The results suggest that WISP1 gene could be important for bone health in humans, as has already been shown in vitro and vivo. The denser fine-mapping, replication, and functional validation will be necessary to understand the mechanisms underlying these associations.

Methods
Subjects. The OOA of Lancaster Pennsylvania are relatively homogenous in terms of both genetic ancestry and lifestyle characteristics. Subjects (n = 1,510) included in this study were adults aged 18 years or older, who participated in the Amish Family Osteoporosis Study (AFOS), the Amish Family Longevity Study (AFLS) and Pharmacogenomics of Anti-Platelet Intervention (PAPI). The protocols and procedures for those studies were approved by the Institutional Review Boards of the University of Maryland and all subjects gave written informed consent. All methods were performed in accordance with the relevant guidelines and regulations. Details of these studies design, recruitment, and phenotyping had been described previously [33][34][35][36] . Briefly, the AFOS was started in March 1997 to study genes that are important for bone health. This was a population-based study to identify individuals with osteoporosis. After the identification of individuals with osteoporosis by BMD, family members were recruited including first-degree relatives aged 20 years and older. Spouses were also invited to participate in the study. The goal of AFLS is to identify genes related to a long and healthy life and to understand the function of these genes. The goal of the PAPI study was to understand the reason why some people do not respond to commonly used medications to prevent myocardial infarction, including aspirin and clopidogrel (Plavix). The subjects in the three studies had bone mineral density (BMD) measured by DXA.

Body composition and bone mineral density. Examinations were conducted at the Amish Research
Clinic in Strasburg, PA 31,33,37,38 . Height was measured by using a stadiometer. Height and weight were recorded without shoes. Body mass index (BMI) was calculated as weight in kilograms divided by height in meters squared. Areal BMD (grams per square centimeter) was measured by DXA at the total hip, hip femoral neck, hip intertrochanter, hip trochanter bone and lumbar spine (L1-L4), using a Hologic 4500 W (Bedford, MA, USA) by a technician certified in bone densitometry. Our study focused on multiple densitometric phenotypes that we considered clinically relevant. For femoral neck, there were 1510 OOA examined (male n = 715 and female n = 795), but for spine, five subjects were excluded, due to either prior spine surgery or structural abnormalities, leaving male n = 713 and female n = 792 participants.
Genotyping and single nucleotide polymorphism (SNP) Selection. Study participants were genotyped using either the Affymetrix 500k or Affymetrix v6.0 SNP chips by the Genomics Core Laboratory at the University of Maryland. SNPs with a minor allele frequency (MAF) ≥1%, a call rate exceeding 95% and conforming to the expectations of Hardy-Weinberg equilibrium (p > 10 −6 ) were used for imputation with IMPUTE-2 using 1000 G CEU reference sample phase2. SNPs with Imputation quality score (INFO) ≥30% were considered. Selection of the tagSNPs was performed based on the OOA genotyping data. Using the aggressive tagger mode of Haploview version 4.2 (http://www.broadinstitute.org/haploview/), we selected 58 tagSNPs which cover all common genetic variation within 5 kb upstream and downstream to WISP1 gene (Chr8: 134198282-134248933, GRCh37.p13) ( Table 6). Association analysis including these 58 tagSNPs was performed.
Statistical Analysis. Summary statistics of baseline clinical characteristics were expressed as unadjusted means ± standard deviations (SD) using the SPSS statistics version 23 (IBM Corporation, N.Y., NY, USA). The association analyses were carried out using in-house software called MMAP (https://mmap.github.io/). The polygenic component was modeled using the relationship matrix derived from the complete 14-generation Amish pedigree structure. We included family structure, study, age, sex, age*sex, as covariates in the association models. BMI was associated with BMD on univariate analysis and was therefore included as a covariate in model 2. Subgroup analyses to determine whether there were differences in gender were performed. Estimation of the additive genetic heritability follows basic quantitative genetic theory, which models the phenotypic covariance (conditional upon covariate effects) between two individuals in a pedigree as a function of their degree of biologic relatedness. Maximum likelihood methods were used to estimate the values of the parameters, such as heritability, that resulted in highest likelihood obtained across all of the pedigrees. Covariates for BMD heritability analysis were study, sex and age. P-values less than 0.05 were considered as significant. Correction for multiple testing was performed using the Bonferroni method for the number of SNPs and traits tested (P = 0.05/ (58 × 5) = 1.72 × 10 −4 ).  Bioinformatics analysis. Analysis of linkage disequilibrium (LD) statistics (r 2 ) surrounding variants of interest was performed using Haploview version 4.2 (http://www.broadinstitute.org/haploview/). Prediction of histone marks and DNAse hypersensitivity sites was performed using HaploReg v4.1 39 , and the five SNPs were annotated in regulatory elements cataloged in Encyclopedia of DNA Elements (ENCODE) project according to UCSC Genome Bioinformatics website (http://genome.ucsc.edu/). The eQTL analyses were performed in GTEx (https://www.gtexportal.org/home/), ExSNP (http://www.exsnp.org/), and primary osteoblasts (obtained from bone biopsies) 21 .