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.

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.

RESULTS

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 improved, 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).

Table 1.

20 representative biological process GO terms associated in the DS trait of 4,659 Korean women (PGS ＜ 0.05)

Set namea	Gene	PGS	Top 3 genes
Extracellular structure organization	150	0.003579	CDH1, WNT7A, NRXN3
Calcium ion transport	127	0.01047	TRPV5, TRPV6, RYR3
Oxygen transport	13	0.008245	HBE1, HBG2, HBD
Telencephalon cell migration	21	0.018368	DAB1, DRD1, LHX6
Glutamate signaling pathway	21	0.026981	APP, GRIK4, GRIN2B
Regulation of neurological system process	142	0.026933	NLGN1, PARK2, GRIN2B
Regulation of synapse structure and activity	22	0.022703	WNT7A, APP, PTK2
Learning or memory	103	0.025144	DBH, AMPH, APP
Regulation of Rho protein signal transduction	95	0.023949	FGD4, TIAM1, FARP1
Dendrite development	29	0.021385	APP, KLHL1, GRIN3A
Startle response	14	0.021831	PARK2, GRIN2B, GRIN3A
Neuron projection morphogenesis	190	0.024503	NRXN3, DMD, CHST3
Sprouting angiogenesis	11	0.026282	BMP4, CDH13, ANGPT1
Central nervous system neuron development	28	0.028691	GNAQ, PTK2, SOX1
Cellular response to heat	16	0.048109	FGF1, MYOF, RBBP7
Filopodium assembly	14	0.04967	FGD4, MTSS1, ARHGAP26
Activation of protein kinase activity	106	0.047758	DGKH, NRG1, BMP4
Response to isoquinoline alkaloid	16	0.046919	GNAO1, SRR, OPRM1
Post-embryonic development	63	0.0487	BMP4, DHCR7, GNAQ
Regulation of GTP catabolic process	116	0.048343	FGD4, AMPH, GNAO1

aThe terms common with the MS results (Table 2) are shown in bold.

Table 2.

30 representative biological process GO terms associated in the MS trait of 4,659 Korean women (PGS ＜ 0.05)

Set namea	Gene	PGS	Top 3 genes
Calcium ion transport	127	7.19E-05	TRPV5, TRPV6, GRIN2A
Neuron projection morphogenesis	190	1.27E-04	CNTN4, ULK2, CNTN2
Cell-substrate adhesion	94	1.24E-04	BVES, ITGA6, ITGB5
Cell morphogenesis involved in neuron differentiation	183	0.002851	CNTN4, ULK2, CNTN2
Adult walking behavior	19	0.002844	SCN8A, SCN1A, CACNB4
Regulation of cell morphogenesis	141	0.002631	VEGFA, GNA12, MAPT
Homophilic cell adhesion	128	0.007898	CELSR1, PTPRT, CDH11
Immune response-regulating cell surface receptor signaling pathway	41	0.009603	SYK, CD247, CACNB4
Activation of immune response	87	0.009152	CFHR1, SYK, CD247
Extracellular structure organization	150	0.013693	NLGN1, TNFRSF11B, ENG
Regulation of cell adhesion	140	0.013691	TESC, NRG1, ENG
T cell receptor signaling pathway	21	0.014998	CD247, CACNB4, FYN
Learning or memory	103	0.015426	HTT, GRIN2A, PARK2
Positive regulation of developmental growth	17	0.015071	MAPT, NTN1, NTRK3
Regulation of positive chemotaxis	21	0.017023	VEGFA, CDH13, FGF10
Filopodium assembly	14	0.020215	ITGA6, ARHGAP26, FGD4
Neurotransmitter transport	80	0.023406	NLGN1, LIN7A, SLC6A11
Regulation of Rho protein signal transduction	95	0.024437	BVES, ITSN1, SYDE2
Microspike assembly	15	0.024645	ITGA6, ARHGAP26, FGD4
Aromatic compound catabolic process	21	0.029948	EPHX1, IDO1, IDO2
Response to alkaloid	54	0.031088	CNR1, SDC2, RYR2
Epithelial cell differentiation	156	0.033825	VEGFA, ENG, KCNE1
Integrin-mediated signaling pathway	56	0.037057	ITGA6, ITGB5, SYK
Startle response	14	0.038299	GRIN2A, PARK2, GRIN2B
Wnt receptor signaling pathway, calcium modulating pathway	20	0.046932	WNT16, WNT9B, FZD2
Epithelial cell proliferation	22	0.048362	TNFSF11, HNF1B, PDGFB
Regulation of cell migration	200	0.047152	VEGFA, GREM1, LAMA4
Regulation of neuron differentiation	156	0.04748	NLGN1, CNTN4, CDK5RAP3
Protein oligomerization	180	0.047725	TRPV5, NLGN1, TEK
Detection of chemical stimulus	33	0.048317	ENG, KCNMB2, RYR2

aThe terms common with the DS results (Table 1) are shown in bold.

The following 10 significant GO terms were replicated between the DS and MS traits: “extracellular structure organization”, “calcium ion transport”, “telencephalon cell migration”, “learning or memory”, “regulation of Rho protein signal transduction”, “startle response”, “neuron projection morphogenesis”, “filopodium assembly”, “response to isoquinoline alkaloid”, and “regulation of GTP catabolic process”. None of them have been reported in previous pathway analyses of BMD GWAS (10, 11).

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-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). Mutations in 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 functional 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 Ca2+ 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 et al. (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.

Fig. 1.

Comparison of -log-transformed P-values of the member genes of the GOs commonly associated between two different skeletalsites: (A) “extracellular structure organization”, and (B) “learning or memory”.

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.

MATERIALS AND METHODS

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 cohorts, 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, m 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.