, INTRODUCTION
Traditionally, measurements of bone mineral density (BMD) are considered essential for diagnosing osteoporosis, which is defined as BMD ≤−2.5 in the T-score [1]. This definition indicates that low BMD significantly increases fracture risk. However, more than half of fragility fractures occur in non-osteoporotic women [2,3]. Thus, other determinants of bone strength, such as the bone microarchitecture and bone geometry, contribute to skeletal fragility [4]. From this perspective, trabecular bone score (TBS), an analytical tool used to measure grey-level textures on lumbar spine (LS) dual X-ray absorptiometry (DXA) images, was introduced to provide information related to the bone trabecular microarchitecture. Previous studies demonstrated that independent of BMD, TBS was associated with the incidence of new fracture and had an additive effect on predicting fracture risk [5,6].
Genome-wide association studies (GWAS) have revealed a large number of loci associated with BMD and fractures. Initial studies identified loci associated with BMD, such as RANKL, OPG, ESR1, FTO, RUNX2, and LRP5 [7-11]. Since then, many GWAS have been performed, revealing more than 100 candidate genes for osteoporosis [12-14]. Among them, novel variants such as ZBTB40, AXIN1, and WLS were selected for further studies [13,15]. However, the genetic determinants affecting trabecular bone microarchitectures and BMD may differ. Some studies observed differences in genetic heritability between cortical and trabecular bones. The heritability of trabecular volumetric BMD has been reported to be as high as 59%, whereas that of cortical volumetric BMD was relatively low at 40% [16]. Additionally, GWAS revealed distinct genetic associations in the lumbar and hip BMD, suggesting that the genetic influences differ in cortical and trabecular bones [13]. However, the genetic contributors to the trabecular bone microarchitecture that influence TBS are unknown.
This study was conducted to identify the genetic determinants of the bone microarchitecture analyzed by TBS. Using both discovery and validation cohorts, we first aimed to identify genome-wide significant genetic variants associated with TBS and then to validate the identified variants in the validation cohort.
METHODS
Participants in the discovery set
Participants were recruited from the Ansung cohort based in Ansung city, Gyeonggi-do, Korea, as a part of the Korean Genome Epidemiology Study (KoGES) [17,18]. The investigation was performed in 2008 and involved 4,206 people in a community-dwelling population. Participants without BMD and TBS (n=200) or GWAS data (n=1,427) were excluded. Men under the age of 50 years (n=168), premenopausal women (n=262), participants being administered anti-osteoporotic drugs (n= 169), or those diagnosed with malignancy (n=15) were also excluded. The total number of participants included in the analysis was 2,451 (1,000 men and 1,451 women). The present study was approved by the Institutional Review Board of Seoul National University Hospital (IRB number 1805-086-946) and was performed according to the Helsinki declaration on the use of human subjects for research. Informed consent was obtained from all patients in the cohort.
Participants in the validation set
Participants in the validation cohort belonged to the Gene-Environment Interaction and Phenotype (GENIE) cohort, which is a sub-cohort of The Health and Prevention Enhancement of the Seoul National University Hospital Gangnam Center in Korea [19]. The cohort includes patients who visited the Seoul National University Hospital Healthcare System Gangnam Center for health check-up between January 2014 and December 2014. Among the 8,000 patients in the GENIE cohort, participants without BMD and TBS data (n=4,085) were excluded. Also, men under the age of 50 years and premenopausal women (n= 1,138), or participants diagnosed with malignancy (n=44) were excluded. In the final analysis, 2,733 (949 men and 1,784 women) were eligible for validation according to the criteria described above for the discovery set. This retrospective study was approved by the Institutional Review Board (H-1712-067-906) of the Seoul National University Hospital; the requirement for informed consent was waived by the institution.
Measurements of BMD and TBS
The BMD (g/cm2) of skeletal sites (LS, femoral neck, and total hip [TH]) and muscle mass were measured using DXA (GE Lunar Prodigy, GE Healthcare, Little Chalfont, UK) and analyzed (Encore Software version 11.0, Encore Software, Minneapolis, MN, USA) according to the manufacturer’s instructions. The BMD precision error (% coefficient of variance [CV]) was 1.7%, 1.8%, and 1.7% for the LS, femoral neck, and TH, respectively, in the Ansung and GENIE cohorts [20,21]. For the LS BMD, the L1-4 value was used for analyses. LSs with compression fracture or severe sclerotic change were omitted from the analyses (e.g., L2-4 were included in analyses if there was a compression fracture at L1). All TBS measurements were retrospectively performed using TBS iNsight Software version 2.0.0.1 (Med-Imaps, Pessac, France) using spine DXA files from the database to ensure that all investigators were blinded to all clinical parameters and outcomes. The software used the raw DXA image of the anteroposterior spine for the same region of interest as the BMD measurement. The precision error (%CV) was 1.5% for the TBS. The instruments were calibrated using anthropomorphic phantoms.
Genotyping and imputation
We extracted genomic DNA from peripheral leukocytes of the participants. The discovery and validation cohorts were analyzed using Affymetrix Genome-Wide Human single nucleotide polymorphism (SNP) Array 5.0 and Affymetrix Axiom KORV1.1-96 Array (Affymetrix, Santa Clara, CA, USA), respectively, which was performed by DNAlink Inc. (Seoul, Korea). Only unrelated subjects with <5% missing genotype were included in the analysis. Markers showing significant deviations from the Hardy-Weinberg equilibrium (P<1.0×10−6), a genotype call rate <0.95, and minor allele frequency <0.01 were excluded. A total of 351,669 SNPs in the discovery cohort and 345,072 SNPs in the validation cohort were directly genotyped and were available for analyses. Imputation was performed using IMPUTE software (https://mathgen.stats.ox.ac.uk/impute) [22]. The HapMap phased genotype information of Japanese individuals from Tokyo, Japan (JPT) and unrelated Han Chinese individuals from Beijing, China (CHB) (build 36 release 22) was used as a reference. Imputed SNPs with a high genotype information content (info >0.5) were used [23]. The same quality control criteria were applied to the imputed SNPs and genotyped SNPs.
Statistical analyses
Association testing was performed using PLINK v1.9 (http://pngu.mgh.harvard.edu/purcell/plink/) [24]. Linear regression assuming an additive genotypic model was used for TBS. Age, sex, and body mass index (BMI) were used as covariates. The results of the discovery cohort can be downloaded from the following webpage: http://bri.snuh.org/bench/_/notice/5250/view.do. In the validation cohort, only variants that passed the significance threshold (P<5.0×10−8) or false discovery rate (FDR) (P<0.05) in the discovery cohort were analyzed [25,26]. The significance threshold for the validation cohort follow-up was P<0.05. Genome-wide significance was considered as P<5.0×10−8. Analyses were performed separately for each sex to identify sex-specific effects.
All analyses were performed with the R language and environment for statistical computing (R Foundation for Statistical Computing, Vienna, Austria). We used Plink v1.9 [24] and data from the 1000 Genome Project Phase 3 [27] to identify proxy SNPs (r2<0.7). Each GWAS locus was defined based on the positions of the left and rightmost proxy SNPs. All SNP and gene locations were relative to the hg38 genome assembly. We also included the closest gene up- and down-stream.
RESULTS
Characteristics of the included cohorts
Anthropometrics and bone phenotypes of the discovery and validation cohorts (KoGES for discovery, GENIE for validation) are shown in Table 1. In the discovery cohort, the mean age of the participants was 62.6 years, and 59.2% of them were female. LS, femur neck [27], and TH BMDs were higher in men than in women. In the validation cohort, as the mean age was 51.1 years and 65.3% of them were women, participants were younger and more likely to be women. The BMI of the subjects was also higher in the validation cohort. For BMD, the LS BMD was higher in the validation cohort than in the discovery cohort (1.184±0.170, 1.007±0.194 g/cm2, P=0.005, respectively). Additionally, in women, FN, TH BMD, and TBS were higher in the validation cohort than in the discovery cohort.
GWAS in the discovery cohort
To detect associations between the genetic variants and TBS, we performed linear regression analysis after adjusting for age, sex, and BMI. In the discovery cohort, 351,669 SNPs were genotyped and passed our quality control filters. After imputation, 1,418,709 SNPs were available for analyses. The Manhattan plot showing significant SNPs is presented in Fig. 1. Two independent variants were selected based on their suggestive associations that passed the significance threshold or FDR post hoc analyses (P<0.05) (Table 2). In the total population, the rs62289585 variant (gene LOC107986169, intronic variant, position chr3: 189493025 [reference GRCh38.p12]) in the intron near trpG1 (anthranilate synthase subunit II) was associated with TBS (P=7.391.0×10−8) and remained significant after adjusting for age, sex, and BMI (P=1.19×10−7). Additionally, the rs55885339 variant (gene CDH4, intronic variant, position chr20:61652860 [reference GRCh38.p12]) was associated with TBS (P=1.65×10−6) and was significant after adjustments (P=3.45×10−6).
Because the TBS is largely influenced by sex, subgroup analysis according to sex was performed. In men, four SNPs were significantly related to TBS, which were rs7206379 (gene WWOX, intronic variant, position chr16:78260046), rs78861306 (gene LOC102723560, intronic variant, position chr16:62751082), rs3851758 (near gene SNX29P2, intronic variant, position chr16: 29248901), and rs1815994 variants (near gene IRX3, intronic variant, position chr16:54523817) in chromosome 16 with reference of GRCh38.p12. In women, five SNPs were also found to be significant, which were rs11630730 (gene MAP2K, intronic variant, position chr15:67565508), rs67059207 (gene SPAG7, intronic variant, position chr17:4965798), rs9507502 (gene MUP58, intronic variant, position chr13:25345858), rs2244721 (near gene KRTAP27, intronic variant, position chr21:30350740), and rs785021 variants (near gene LOC101928442, intronic variant, position chr15:47856460) with reference of GRCh38.p12. Significant SNPs in men and women differed from those found in the total participants.
GWAS in the validation cohort
Eleven variants were further evaluated in the GENIE cohort, which was used as the validation cohort (Table 2). Among them, two SNPs showed significant associations with TBS in the validation cohort (P<0.05) (Table 2). In men, the rs1815994 variant in the intron of IRX3 was associated with TBS (P=3.74×10−5 in the discovery cohort, P=0.027 in the validation cohort). In women, the rs11630730 variant in the intron of mitogen-activated protein kinase kinase 5 (MAP2K5) was related to TBS in both cohorts (P=3.74×10−5 in the discovery cohort, P=0.027 in the validation cohort). Regional association plots of the variants near rs1815994 and rs11630730 are shown in Fig. 2.
Phenotypes related to bone according to the rs1815994 and rs11630730 variants
Among the two variants which passed the threshold for the discovery and validation cohorts, men with the AA genotype of rs1815994 variant near IRX3 showed a higher TBS than those with the AG or GG genotype (1.438±0.003, 1.414±0.010, 1.359±0.030; P=0.038, P=0.010, respectively) (Supplemental Fig. S1). There was a similar trend in the LS BMD, in which men with the AA genotype showed a higher LS BMD than those with the AG or GG genotype (1.101±0.009, 1.050±0.019, 0.962±0.097 g/cm2; P=0.041, P=0.020, respectively). However, the femur neck and TH BMD were similar among genotypes of the rs1815994 variant.
According to the rs11630730 variant near MAP2K5, women with the AA genotype showed a higher TBS than those with the AG genotype (1.394±0.002, 1.346±0.017, respectively; P= 0.037) (Supplemental Fig. S2). The LS BMD was higher in women with the AG genotype than in those with the GG genotype (0.953±0.007, 0.707±0.090 g/cm2, respectively; P=0.046). The femur neck and TH BMDs were similar among genotypes of the rs11630730 variant. The detrimental effects of the rs1815994 variant in men and rs11630730 variant in women were also identified in association analysis (β=−0.0281, β=−0.0465, respectively) (Table 2).
In the validation cohort, men with the AG genotype of rs1815994 variant had higher TBS than those with GG genotype (1.449±0.066, 1.399±0.065, P=0.042, respectively). Also, men with AA and AG genotype of rs1815994 variant higher LS BMD than those with GG genotype (1.249±0.184, 1.261±0.168, 1.071±0.161, P=0.032, respectively). Women with AG genotype of rs11630730 variant had higher TBS than those with GG genotype (1.391±0.065, 1.312±0.087, P=0.039, respectively). The femur neck and TH BMD were similar between genotypes in both genders.
DISCUSSION
In this 2-stage GWAS, we discovered and confirmed that rs1815994 near IRX3 in men and rs11630730 near MAP2K5 in women were associated with the deterioration of bone microarchitecture reflected by TBS in both cohorts. We also identified additional novel variants such as rs7206379, rs78861306, and rs3851758 near WWOX, LOC102723560, and SNX29P2, respectively, in men and rs67059207, rs9507502, rs2244721, and rs785021 variants near SPAG7, MUP58, KRTAP27, and LOC101928442, respectively, in women. However, these genes were significant only in the discovery cohort. This is the first study to identify a presumable functional variant associated with TBS.
IRX3 encodes for iroquois-class homeodomain protein IRX-3, which is a member of the iroquois homeobox family and a transcriptional factor involved in neural development [28]. IRX3 is reported to be related to obesity [29], angiogenesis [30], and limb skeletal formation [31]. However, the functional role of IRX3 in the bone microarchitecture is not fully understood. A previous study showed that Irx3 and Irx5 regulate early skeletal formation by modulating sonic hedgehog signaling [31]. As sonic hedgehog signaling is known to alter the architecture and mechanical properties of trabecular bone [32], the rs1815994 variant near IRX3 may regulate the bone microarchitecture in older adults. In another study, Irx3 and Irx5 were found to be associated with skeletal mineralization [33]. Osteoblast-specific deletion of Irx3 and Irx5 causes mineralization defects in cranial bones, as well as bone fragility similar to Hamamy syndrome, which is a rare autosomal recessive syndrome involving mental retardation and osteopenia with repeated fractures [34]. Additionally, a recent study showed that Irx3 regulated chondrogenic differentiation in mouse mesenchymal cells. Irx3 expression was increased with chondrocyte differentiation of mesenchymal cells following Bmp2 treatment in the study [35]. As reported, chondrocyte-induced signaling such as Indian hedgehog was essential for maintaining the growth plate and articular surface and was required for sustaining the trabecular bone structure [36]. Therefore, IRX3 near rs1815994 variant may be a novel target for maintaining trabecular microarchitecture, which should be further examined.
Another notable result of this study was the identification of the rs11630730 variant near MAP2K5 associated with TBS in women. MAP2K5 encodes for MAPK/ERK kinase 5 and is ubiquitously expressed in most human tissues [31]. This gene plays a critical role in activating the MAPK7/ERK5 pathway, which is essential for skeletal development and homeostasis [37]. In a previous study, the MAPK/ERK pathway was highly activated in osteoblasts in mature osteoblasts in vitro and in vivo. When the ERK pathway was inactivated by knock-down of MAP2Ks, mice without MAP2K activity showed severe osteopenia with growth retardation. Interestingly, inactivation of the ERK pathway in adult mice with tamoxifen-inducible Cre recombinase reduced the trabecular bone mass, number, and thickness in mice [37]. In this regard, modulating MAP2K5 activity by the rs11630730 variant may have a significant role in maintaining the trabecular bone microarchitecture. Other studies showed that MAP2K5 is essential for the early stages of myogenesis [38-40]. As ERK5 and MAP2K5 modulate IGF-2 to regulate promyogenic activities, muscle atrophy by an impaired MAPK5/ERK5 pathway can cause impaired bone microarchitecture [39,41]. However, further studies are needed to confirm the association of rs11630730 with decreased expression of MAP2K5 in the bone or muscle. As patients with different genotypes had different BMD and TBS, specific genotypes of the variants could be used as one of the indicators for patients with a low risk of degraded microarchitecture or osteoporosis. In specific, the AA genotype of rs1815994 variant in men and the AA genotype of rs11630730 variant in women could be potential protective markers for osteoporosis or conserved microarchitecture.
Other genes may have clinical implications as markers for the bone microarchitecture, although they were significant only in the discovery cohort. The rs55885339 variant near CDH4 showed a significant association in the total population in the discovery cohort but was not significant in the validation cohort. CDH4 encodes R-cadherin/cadherin-4 protein, which is a cell adhesion molecule. Cadherins are well-known to be anchored to the cytoskeleton via binding to β-catenin [42]. As β-catenin is vital in the canonical Wnt signaling system, cadherins can modulate osteogenic differentiation, skeletal development, and maintaining homeostasis by regulating binding affinities with β-catenin [43]. In men, we also found an association of the rs7206379 variant, which is near WWOX. WWOX encodes WW domain-containing oxidoreductase, which is a tumor suppressor gene in osteosarcoma and is also related to bone metabolic disease such as osteoporosis [44]. WWOX associates with RUNX2 to suppress its transcriptional activity in osteoblasts [45]. In this regard, both rs55885339 and rs7206379 variants near CDH4 and WWOX may be novel modulators of the bone microarchitecture, which should be validated in further replication studies. Variants and genes found in this study were novel and have not been previously described in other GWAS of BMD or fracture. Additionally, the genetic control of trabecular bone microarchitecture may differ from that of BMD. However, since TBS is a calculated score from BMD results, which is an indirect measurement of the microarchitecture, the result should be carefully interpreted.
Although the variants in this study have not been reported in the previous studies regarding BMD, known variants associated with BMD may be relevant to the variants from the study. In previous reports, over 100 loci, including FTO and RUNX2, were found to be linked with BMD [7-11]. Interestingly, SNPs of FTO which is a known variant associated with BMD and IRX3 from our study were in strong linkage disequilibrium and associated with the risk of obesity [12]. As obesity is correlated to the trabecular microarchitecture of bone [13,14], polymorphisms of IRX3, along with FTO, may have an association with the bone microarchitecture. Also, mitogen-activated protein kinases (MAPKs), including MAP2K5 from our study, are traditional signal transducers and reported to act in phosphorylating RUNX2 [15]. Therefore, variants in MAP2K5 could result in different activities of RUNX2, which is also a well-known variant associated with osteoporosis [16], which may have a combined effect on the bone trabecular microarchitecture. However, future works with a focus on the combined effect of different variants are needed.
We believe that our study provides novel evidence that previous GWAS could not suggest from GWAS with total BMD. This is the first study to demonstrate an association of the variants with TBS, and thus further replication studies are required to confirm and to identify additional variants involved in the trabecular bone microarchitecture. These findings may result in the identification of novel osteoporosis drug targets, particularly for IRX3 and MAP2K5. There were some limitations to this study. A major limitation is the relatively modest sample size for GWAS, which may have caused us to miss some true-positive associations. Additionally, the interpretation of the results needs caution since TBS only reflects the bone microarchitecture, does not measure it directly. However, TBS has been shown to be associated with fracture independently of BMD, and has been widely validated and used in previous clinical studies [46,47]. Importantly, patients in the discovery cohort were older and had a higher BMI than those in the validation cohort. To minimize these differences between cohorts, we adjusted for age, sex, and BMI when calculating the linear association of variants with TBS in both cohorts. However, the different baseline characteristics between cohorts should be considered when interpreting the results, mainly when the variants showed a significant association in only one cohort.
In this GWAS, we identified an intronic variant of IRX3, rs1815994, in men and MAP2K, rs11630730, in women associated with the trabecular bone microarchitecture analyzed based on the TBS. These results broaden the understanding of the genetic factors contributing to the trabecular microarchitecture and may provide novel targets for fracture prevention.
PDF Links
PubReader
ePub Link
XML Download
Full text via DOI
Download Citation
CrossRef TDM
Supplement1
Print
