Gene-set based genome-wide association analysis for the speed of sound in two skeletal sites of Korean women

The speed of sound (SOS) value is an indicator of bone mineral density (BMD). Previous genome-wide association (GWA) studies have identified a number of genes, whose variations may affect BMD levels. However, their biological implications have been elusive. We re-analyzed the GWA study dataset for the SOS values in skeletal sites of 4,659 Korean women, using a gene-set analysis software, GSA-SNP. We identified 10 common representative GO terms, and 17 candidate genes between these two traits (PGS ＜ 0.05). Implication of these GO terms and genes in the bone mechanism is well supported by the literature survey. Interestingly, the significance levels of some member genes were inversely related, in several gene-sets that were shared between two skeletal sites. This implies that biological process, rather than SNP or gene, is the substantial unit of genetic association for SOS in bone. In conclusion, our findings may provide new insights into the biological mechanisms for BMD. [BMB Reports 2014; 47(6): 348-353]


INTRODUCTION
Osteoporosis, which leads to an increased risk of fracture, is a common disease in elderly women (1,2). The tests for osteoporosis or risk of fracture consider multiple factors (3). One of them is bone mineral density (BMD), which is a polygenic trait, with high heritability. Twin and family studies have shown that between 50 and 85% of the variance in peak BMD is genetically determined (4). Until now, 70 genes that may have influence on BMD have been identified, through 17 genome-wide association (GWA) studies (http://www.genome. gov/gwastudies/). However, biological connections of these candidate genes to BMD have been difficult, and little overlap has been observed between the results of those GWA studies. Genes do not work independently; instead, complicated molecular networks and cellular pathways are often involved in disease susceptibility and progression (5,6). Moreover, pathways, rather than individual genes, have high effect (7). Therefore, several gene-set analysis (GSA) methods and tools based on Gene Ontology (GO) or pathway have been developed, in order to facilitate more powerful interpretation for the GWA studies dataset (8). By analyzing SNPs on the gene-set level, GSA was able to find many coordinated association patterns that might be lost by individual marker analysis (9). Recently, pathway analysis for GWA studies has been used to investigate many common diseases or quantitative traits, including breast cancer, Alzheimer's disease, BMD, hypertension, and diabetes mellitus (8). Pathway analyses for BMD identified several significant pathways, such as regulationof-autophagy pathway, and gamma-hexachlorocyclohexane degradation (10,11). While these analyses were performed using GWA results that examined male and female subjects together, it would be desirable to perform GSA on GWA results based on female data, considering the gender difference (2). The speed of sound (SOS) value measured by quantitative ultrasound is an indirect indicator of BMD (12). Here we applied GSA to the GWA results that were performed, using SOS values measured at both Distal Radius (DS), and Mid-shaft Tibia (MS), in Korean women. As the phenotype correlation between these two traits was 0.455 ( Supplementary Fig. 1), the two datasets were analyzed separately, and the results were compared, to discover common biological pathways affecting independent skeletal sites.

Correlation of SNPs, genes, and biological process GO terms between DS and MS traits
We performed GWA analysis for SOS in 4,659 Korean women, and then used the P-values of the imputed 1,827,004 SNPs from the result. All P-values without a cutoff were compared pairwise at the SNP marker level, using scatter plots. The result showed poor correlation, as the correlation coefficient was 0.087 (Supplementary Fig. 2A). The SNP P-values were summarized gene-by-gene, and the best P-value was assigned for each gene. The scatter plots involving all genes showed imhttp://bmbreports.org BMB Reports proved, but still low correlation between two traits (DS and MS), as the correlation coefficient was 0.33 ( Supplementary  Fig. 2B). The gene-level P-values were summarized gene-set by gene-set, using GSA-SNP (9). Supplementary Fig. 2C shows a scatter plot of the log-transformed PGS between DS and MS traits, demonstrating a fairly increased correlation (0.53). The concordance at the gene-set level was somewhat higher, than that at the SNP level and gene level.

Results of GSA using GSA-SNP tool
We summarize the result of gene-set analyses of the GWA results of 4,659 Korean women, for SOS values in two different skeletal sites. In our gene-set analyses, we used biological process terms of the GO database. As these terms are hierarchically arranged, the resulting list of terms is highly redundant. The lists were pruned manually, resulting in 20 and 30 representative biological process GO terms for the two traits, respectively (PGS ＜ 0.05) ( Tables 1 and 2).
One of the genes showing repeatedly dominant association in member genes of several significant GO is the transient receptor potential family, vanilloid type-5 (TRPV5) gene. TRPV5 and TRPV6 enable bone formation, and relate to osteoporosis, by regulating calcium uptake (13).

Genetic factors with significant association
Not all member genes of the gene-sets that are identified from the GSA show significant association with the traits. We focused on the members that were significantly associated (PGWAS ＜ 1 × 10 -3 ), and highly effective (|β| ＞ 200), where β represents the regression coefficient. There were 6 and 11 such genes for DS and MS, respectively (Supplementary Tables  1 and 2); among them only one gene (FGD4) was common. Almost half of these findings have been reported for association with BMD by previous studies (14)(15)(16)(17)(18)(19)(20)(21)(22)(23). For example, ALCAM may play a role in the progress of osteogenic differentiation, because ALCAM+ cells can differentiate into osteoblasts, which can support osteoclastogenesis (14). The Ca2+-like effects of CD38 might be relevant physiologically in the metabolic control of bone resorption via NAD(+). Therefore, CD38/Ca2+/IL-6 pathway may have a critical role in coupling an osteoclast's metabolic activity with its resorptive function (15). Human osteoblasts express a repertoire of cadherins, including N-cadherin, cadherin-11, and cadherin-4 (CDH4). Expression of CDH4 mRNA and protein was strongly induced by dexamethasone in osteoprogenitor marrow stromal cells, and was stimulated in normal human trabecular bone osteoblasts (16). MAPKs function to regulate the key transcriptional mediators of osteoblast differentiation, with ERK and p38 MAPKs phosphorylating RUNX2, the master regulator of osteoblast differentiation. In addition, MAPKs in osteoblasts play a role in the maintenance of bone mass (17  P2RX7 gene have been associated with low lumbar spine BMD, and accelerated bone loss in post-menopausal women (18). FGFR signaling is known to be important in the initiation and regulation of osteogenesis, and the IIIc alternative of Fgfr2 is a positive regulator of bone formation, affecting mainly the osteoblast (19,20). SYK(-/-) progenitors are similarly defective in osteoclast development and bone resorption (21). Expression and activity of TRPM7 gene to calcium ion transport is modulated by extracellular Mg2+ and Ca2+ availability, indicating that TRPM7 channels are involved in intracellular ion homeostasis, and proliferation of osteoblasts (22). Moreover, TRPM7 knockdown inhibits osteoblast proliferation, in mature bone cells (23).

DISCUSSION
The genetic factors for complex human traits discovered by GWAS are usually based on the analyses of individual SNPs. However, GWAS accounts for only a small proportion of the heritability of complex traits, and it is not well suited to the detection of small effects of multiple SNPs (8). On the other hand, GSA reveals the cumulative contribution of the func-tional genes, and is useful in biological interpretation.
In our study, we found 17 candidate genes all together, and 20∼30 representative biological process GO terms that may influence the SOS value of Korean women. The association of the gene-sets was well replicated; whereas, the phenotypes and genetic factors were significantly different between the DS and MS traits. Therefore, GO terms commonly implicated between DS and MS may play important roles in the bone mechanism. For example, "extracellular structure organization" was identified to be the top GO associated with DS in our study. Extracellular matrix (ECM) proteins induce the osteoblast differentiation of human mesenchymal stem cells (24). Osteoblasts create the nano-composite structure of bone, by secreting a collagenous ECM, on which apatite crystals subsequently form (25). "Calcium ion transport" was well known to affect BMD, through other studies. The osteoclast has several other calcium transport proteins, including a Ca 2+ ATPase, which is expressed highly. Calcium signaling promotes osteoclastogenesis, and bulk calcium transport by the osteoblast. Moreover, disorders in calcium removal from bone matrix by osteoclasts cause osteopetrosis (26). In the "regulation of Rho protein signal transduction" and "regulation of GTP catabolic process", Chellaiah http://bmbreports.org BMB Reports  (2000) demonstrated that RHOA protein is essential for osteoclast motility and bone resorption by transducing active and inactive RHOA into avian osteoclasts (27). In addition, the RhoGTPase-RhoGEF pathway plays an important role in bone cell biology and osteoporosis (28). Several GO terms associated with nerve were also identified, as follows: "telencephalon cell migration", "learning or memory", "startle response", and "neuron projection morphogenesis". Nerve fibers with active expression of various neural transmission ligands were demonstrated to be in close spatial association with bone cells. Moreover, receptors for these neural ligands are expressed by bone cells, and administration of these neural transmission molecules has potent effects on bone cells. Metabolic control of bone is influenced by the nervous system, and potential transmitters of this influence include glutamate, calcitonin gene-related protein (CGRP), leptin, and so on. Disorders of nerves -central or peripheral -can have substantial effect on bone health and repair (29). For the "learning or memory", Loskutova et al. (2009) reported that BMD is reduced in the earliest clinical stages of Alzheimer's disease (AD), and associated with brain atrophy and memory decline, suggesting that central mechanisms may contribute to bone loss in early AD (30). Zhang et al. (2001) reported that BMD in the elderly is associated with verbal memory impairment, and that the mechanisms underlying this relation are not understood, but cumulative exposure to estrogen may play a role (31). In the "filopodium assembly", an actin superstructure that links two precursor cells appears temporarily during the secondary fusion of osteoclasts. In one type of secondary fusion, the osteoclasts use a filopodium-like protrusion that linked the cells (32). For the "response to isoquinoline alkaloid", preventing the differentiation and function of osteoclasts at the early stage was an important anti-bone destruction mechanism of Norisoboldine (NOR), as the major isoquinoline alkaloid in Radix Linderae, which might be attributed to the inhibition of ubiquitination of TRAF6, the accumulation of TRAF6-TAK1 complexes, and the activation of MAPKs/NF-κB/c-Fos/NFATc1 pathways (33). We also discovered that member genes showing strong association signals are different, by comparing the P values of the member genes of a common gene-set, between different skeletal sites. See Fig. 1 for representative scatter plots of the gene-sets commonly showing association in two skeletal sites. The "extracellular structure organization" GO, as an example, was significantly associated in both the DS (PGS = 0.00358) and the MS (PGS = 0.01369) traits. As shown in Fig. 1A, a subset of the member genes that included CDH1, an epithelial cadherin from the cadherin superfamily, showed strong signals in the DS trait (P value: 2.05E-07), but not in the MS trait (P value: 0.7784). However, a different subset of member genes that included NLGN1, a member of a family of neuronal cell surface proteins, showed a relatively stronger signal in the MS trait, than that in the DS trait. The association between CDH1 and osteoblast function has been well established in the literature (34). This suggests that pathways, rather than SNP or gene, are the substantial units of genetic association for quantitative traits, such as BMD.
In conclusion, GSA to the GWA studies datasets indicated significant biological processes that may regulate BMD. These findings may be supplementary evidence, and provide new insights into the bone mechanisms.

Study samples and association analysis
The samples and genotype data used in this study have been previously described (35). Briefly, through the Korea Association Resource (KARE) project, a total of 10,038 participants were recruited from Ansan and Ansung population-based cohttp://bmbreports.org horts, aged 40 to 69. 10,004 samples were genotyped, using the Affymetrix Genome-Wide Human SNP array 5.0. A total of 352,228 markers in 8,842 individuals were obtained, after removing samples and markers that failed a quality control test (35). We studied 4,659 women samples of them, and used the P values of the imputed genotypes (1,827,004 SNPs). SNP imputation has also been described (35). Briefly, using the IMPUTE program (36) the KARE genotypes were supplemented, by imputing SNP genotypes based on 90 individuals from the unrelated Chinese in Beijing (CHB), and Japanese in Tokyo (JPT) founders in HapMap. SOS was tested for association by linear regression analysis with dominant, additive, and recessive models, after adjusting for age and height as covariates, using PLINK.

GSA-SNP software
GSA was applied to the SOS GWA studies dataset of Korean women, using the GO database, which is composed of the gene-sets having 10-200 members (2,476 biological process GO terms). The Z-statistic method is applied, by using GSA-SNP (9) with the default options. SNPs residing inside or within 20 kb of the boundary of each gene are gathered, and the second best P value is assigned to the gene. The gene score is calculated as -log of the P value assigned to the gene. The Z-statistic is then performed for each gene-set (GS), where,   is the average of gene scores [-log(k-th beset P)] in a geneset, m0 and σ are the mean and the standard deviation of all the gene scores, and n is the number of genes in the gene-set.
The P values of each gene-set are computed under the assumption of a normal distribution of the Z-statistics, followed by multiple testing correction, using the false discovery rate method.