Integrative genomic analysis for the functional roles of ITPKC in bone mineral density

Osteoporosis is a skeletal disorder characterized by low bone mineral density (BMD) that increases the risk of fracture [1]. Its prevalence among the elderly population in Taiwan is 24.15% in females and 11.03% in males [2]. Several factors contribute to the pathogenesis of osteoporosis, including age, sex, and genetics [3]. The heritability of BMD has been estimated at 50–85% according to twin and family studies [4,5], and numerous genes have been identified as being related to osteoporosis and BMD in candidate gene-based studies and genome-wide association studies (GWASs) [6]. The largest meta-analysis of GWASs, published by the GEFOS consortium, revealed several loci associated with BMD [7,8]. Most of the BMD-related common variants function within the RANK–RANKL–OPG signaling pathway, and the low-frequency variants of EN1 have effects on BMD as well. Based on these and other findings, BMD is widely considered to be a polygenic trait.

Bone loss is due to imbalances between osteoclastic bone resorption and osteoblastic bone production. Store-operated calcium (SOC) entry is an important regulator of osteoclast differentiation [9,10], and when early-stage osteoclasts receive the signals from the receptor activator of the nuclear factor κ-B ligand (RANKL), the store-operated calcium channel will be activated, causing calcium influx [9]. This influx leads to calcium oscillations that activate downstream nuclear factor of activated T cells c1 (NFATc1) to promote osteoclastogenesis [9,11]. The SOC channel is the main calcium influx pathway in non-excitable cells—including T cells, mast cells, and osteoclasts—and it is activated by depletion of calcium stores in the endoplasmic reticulum (ER) [12–14]. Upon depletion of stored calcium in the ER, stromal interacting molecule 1 (STIM1) will aggregate and bind ORAI1 to induce influx of extracellular calcium [15]. Inositol-trisphosphate 3-kinase C (ITPKC) is a negative regulator of the SOC entry, which affects the NFATc1 signaling pathway through the phosphorylation of inositol 1,4,5-trisphosphate (IP₃) [16,17]. Once IP₃ is phosphorylated to IP₄, IP₃ receptors on the ER are incapable of activating SOC entry, contributing to calcium reductions in the cell [18].

In 2003, Mentaverri et al. [19] reported that 2-aminoethoxy-diphenyl borate (2-APB) and SKF-96365, two SOC channel blockers, can significantly decrease bone resorption and the survival of osteoclasts. Additionally, knockdown of STIM1 expression or inhibition of IP₃R in bone marrow macrophages reduced calcium signaling and diminished osteoclast differentiation [20]. Hwang et al. [21] demonstrated that by using short hairpin (sh)RNA to silence ORAI1, a SOC entry component in the plasma membrane, RANKL-induced osteoclastogenesis was blocked. Moreover, compared with wild-type mice, ORAI1^−/− mice lacked multinucleated osteoclasts and exhibited markedly decreased cortical ossification and BMD [22,23]. Together, these reports strongly suggest that SOC entry is an important regulatory signal for bone remodeling.

Since the RANKL-induced calcium influx in the early stages of osteoclastogenesis occurs via SOC entry [24], genes in the SOC pathway are likely to be involved in the regulation of bone metabolism. As such, ITPKC, an upstream regulator of SOC entry, is a potential candidate gene for osteoclastogenesis as well. However, the genetic relationships between STIM1, ORAI1, ITPKC, and BMD remained unclear in humans. In the present study, we aimed to investigate the association between these genes and BMD by examining data from the Gene Expression Omnibus (GEO) database and individual’s genotype from 1044 subjects.

Gene expression datasets were collected by searching the GEO database (https://www.ncbi.nlm.nih.gov/geo/) with the following key words: ‘osteoporosis’, ‘BMD’, and ‘Homo sapiens’. We included studies in the analysis of their phenotype data that are available. Other diseases as their major phenotype or the study only contained control samples were excluded. For each study, we extracted information including sample number, platform, phenotype, cell type, and gene expression data.

Raw expression values of each dataset were normalized using Robust Multi-array Average (RMA) algorithm [25]. If raw CEL files were not available, the processed expression values were downloaded directly form GEO. Bioconductor hgu133a.db package was used to annotate each probe ID with its corresponding gene symbol. Each dataset was normalized and mapped to gene symbol individually, and we conducted a meta-analysis by combing expression mean and standard deviation value of each study afterward. R meta package was used to perform meta-analysis and generate forest plot. A statistical test of heterogeneity between studies was estimated by I². Both the fixed-effect and random-effects model were implemented in this analysis. We used meta-regression under REML model to analyze the high heterogeneity result. The metafor package of R was used to conduct this analysis.

We enrolled a total of 1044 patients from Wan-Fang Hospital, Taipei, Taiwan. The study protocol conformed to the Declaration of Helsinki. Males and postmenopausal women aged ≥55 years who visited Neurological Clinics of Wan Fang Hospital due to back pain or lumbar radiculopathy were recruited for the present study. Patients with pathological fractures or high-impact fractures (such as those due to motor vehicle accidents) and continuous steroid use (of over 6 months) were excluded. Patients with the long-term inflammatory disease were also excluded.

BMD was measured by dual energy radiograph absorptiometry with standard protocols at the lumbar spine (LS; L2-4 or L1-4) and femoral neck (FN). Vertebral fractures were assessed by digital measurements of morphologic changes on a lateral radiograph of the thoracolumbar spine. We collected clinical information of the subjects such as age, gender, and BMI. The study was approved by the Joint Institutional Review Board of Taipei Medical University. All subjects provided written informed consent.

DNA was extracted from whole-blood samples using a Gentra (Qiagen, Valencia, CA) extraction kit and 70% alcohol precipitation. The region we used to select tagging SNPs was defined by the position of each gene locus ± 1500 KB. Based on HapMap SNP database (release27 phase II + III Feb 09, dbSNP b126) and the Haploview 4.2, our study selected tagging SNPs among each linkage disequilibrium (LD) block (r² = 0.8). Besides, the SNP with any prior functional report had the priority to be selected. We selected six tagging single-nucleotide polymorphisms (tSNPs) of ITPKC, four tSNPs of STIM1, and five tSNPs of ORAI1 with a minimum allele frequency of ≥10% in a Beijing Han Chinese and Taiwanese population (Supplementary Figures S1–3). Genotyping was performed with a TaqMan Allelic Discrimination assay (Applied Biosystems, Foster City, CA). A polymerase chain reaction (PCR) was carried out using an ABI StepOnePlus Thermal Cycler (Applied Biosystems). In a subsequent PCR, the fluorescence from specific probes was detected and analyzed through the System SDS software version 2.2.2 (Applied Biosystems).

To investigate the association between gene expression profiles and the SNPs of ITPKC, we also queried the GTEx Portal (http://www.gtexportal.org/home/) that includes a variety of tissue expression quantitative trait loci (eQTLs). HaploReg browser (www.broadinstitute.org/mammals/haploreg), which provides regulatory elements estimated by ENCODE and Epigenomics project data, was also used to discover the potential influence of these SNPs.

R 3.2.0 was used for the statistical analyses. Associations between genotypes and BMD at the two sites were tested by the likelihood ratio test. Age, body mass index (BMI), and sex were adjusted in the models as potential confounders. Pairwise LD among genotyped SNPs was assessed and used to define haplotype blocks via Haploview software vers. 4.1 [26]. The general linear model (GLM) implemented in R were used to examine the associations between haplotypes and BMD. To correct for multiple testing, the false discovery rate (FDR) was applied, and q values were estimated to control for proper type I errors. FDR q values of <0.05 were considered statistically significant.

After carefully curating available data and removing duplicate datasets, five distinct gene expression datasets were downloaded from GEO. Three datasets were based on RNA extracted from peripheral blood monocytes, and the other two datasets were derived from RNA collected from B lymphocytes and mesenchymal stem cells, respectively. Overall, we included 155 individuals from five datasets. Among them, 79 belonged to the high BMD group and 76 belonged to the low BMD group. Most of the subjects were female, except one male in the high BMD group. The detailed parameters of each study are displayed in Table 1.

Figure 1 shows the forest plot of ITPKC expression profiling for the low and high BMD groups across all examined studies. The mean difference in ITPKC expression between low BMD and high BMD groups was −0.07 ( 95% CI: −0.12 to −0.01). Based on a fixed effect model, the expression level of ITPKC was significantly associated with BMD (p = 0.03), even under a random-effects model (p = 0.04). There was no severe heterogeneity across studies (I² = 40%). Neither STIM1 nor ORAI1 showed any significant association between expression level and BMD status (Supplementary Figures S4 and S5). The heterogeneity of STIM1 was 89%, so we further conducted a meta-regression. This high value was mainly due to differences in cell types, which accounted for 94.93% of the heterogeneity.

A total of 1044 individuals (794 females and 250 males) were recruited in this study (Table 2). The mean age ± standard deviation was 68.5 ± 9.4 years for females and 71.2 ± 9.8 years for males.

The frequencies of the tag SNPs for ITPKC (rs7257602, rs7251246, rs890934, rs10420685, rs2607420, and rs2290692), STIM1 (rs2304891, rs3750996, rs1561876, and rs3750994), and ORAI1 (rs12320939, rs12313273, rs7135617, rs6486795, and rs712853) were similar to those from the Taiwan Biobank [27] (Supplementary Table S1). After adjusting for age, sex and BMI, the SNP, rs2607420, was significantly associated with LS BMD (q = 0.028). Furthermore, SNP rs10420685 was associated with FN BMD (p = 0.028) (Table 3); however, the association did not reach statistical significance after multiple testing corrections (q = 0.131). The genetic polymorphisms in STIM1 and ORAI1 did not reach statistical significance in either sites (Supplementary Table S2).

In order to elucidate the most important haplotype of ITPKC, we further constructed a LD map for ITPKC. Two haplotype blocks were found after pair-wise LD analysis of the genotyped SNPs (Figure 2). A haplotype association analysis was then performed on each block. For haplotype block one, formed by rs7251246, rs890934 and rs10420685, the C-G-A haplotype was significantly associated with LS BMD compared with the reference T-T-A haplotype after adjusting for covariates (p = 0.008) (Table 4). The C-G haplotype of block two was compared with the T-C haplotype and significant association was found with LS BMD (p = 0.008). Neither haplotype block showed a significant association with FN BMD (Table 4).

In the present study, we first conducted a meta-analysis of STIM1, ORAI1, and ITPKC gene expression based on GEO data. We then further investigated the association between BMD and SNPs in these three genes using 1044 subjects. Our results revealed that among the genes we examined, only ITPKC expression was positively correlated with BMD. In addition, one intronic SNPs of ITPKC, rs2607420, showed strongly association with LS BMD, yet there was no statistically significant association between polymorphisms of ITPKC and FN BMD (Table 4).

Previous studies have illustrated the importance of ORAI1 in bone metabolism. In the present study, we found that the genetic association between ITPKC and BMD is more pronounced. This observation may due to the fact that ITPKC acts as a negative regulator of the SOC channel, which can directly influence calcium signals. Therefore, our data suggest that either the expression level or genetic polymorphisms of ITPKC may be a better biomarker to predict BMD compared with ORAI1 or STIM1-related measurements.

Our results also showed that SNPs may have site-specificity effects, which is consistent with previous studies [7]. Because the proportions of constituent cortical and trabecular bones are distinct at different body locations [28], and various genes differentially regulate the two types of bones, the same genetic polymorphism may have different impacts at different sites [29,30]. Analysis of eQTLs using the GTEx Portal indicated that rs2607420 was significantly associated with the expression of ITPKC, ADCK4, and C19orf54 (Supplementary Table S3). According to 1000 Genomes (Pilot 1 CHB + JPT), rs2607420 is in high LD (r² > 0.6) with other SNPs that are located in C19orf54 and ADCK4 (Supplementary Figure S6). Moreover, the HaploReg browser indicated that rs2607420 is an intronic SNP, which lies in enhancer histone marks (H3K27ac and H3K9ac) and DNase I-hypersensitivity site in several cell lines including primary T cells and B cells from peripheral blood [31]. These results support a functional role for ITPKC in BMD regulation.

A limitation of our study is that we did not collect bone remodeling markers and fracture data to conduct association analyses. These types of data may be required to investigate the further relationships between genetic polymorphisms and osteoporosis-related clinical events. Besides, by using tagging SNP study design we did not genotype all the genetic variants of ITPKC. Target sequencing may be used to detect the potential causal variants. Importantly, other replication studies are needed to validate our results. Although several bioinformatics tools have allowed us to understand the potential functions of ITPKC polymorphisms in determining the risk of osteoporosis, functional studies on animal models in relevant tissues are still needed to clarify the underlying mechanisms.

In summary, our results not only revealed a relationship between ITPKC expression and BMD, but also identified specific BMD-related loci within the ITPKC gene. These results highlight the role of ITPKC in determining inter-individual variations in BMD. Our findings may be useful in the development of novel diagnostic tools or treatment targets for osteoporosis in the future.

This work was supported by grants from Ministry of Health and Welfare [MOHW103-TDU-B-212-113001 and MOST 105-2628-B-038-001-MY4]; Wan Fang Hospital [105-wf-eva-13]; the National Science Council, Taiwan, ROC [MOST 104-2320-B-038-016 and MOST 105-2314-B-038-028]; the Higher Education Sprout Project by the Ministry of Education in Taiwan [MOE DP2-107-21121-A-0]; and Taipei Medical University [12310-0223].

The authors declare that there are no competing interests associated with the manuscript.

Experimental design: W.C.C. and H.F.L.; Implementation of experiments: H.F.L. and Y.W.H.; Data analysis: H.F.L. and H.S.C.W.; Contribution of reagents/materials/analysis tools: K.S.H., B.K.C., E.C.C. and W.C.C.; data interpretation: W.C.C., K.S.H., H.T.L., H.F.L., and S.I.; Clinical samples, data collection: K.S.H., H.T.L.; Revision of manuscript: H.F.L., S.I, E.C.C., and W.C.C.; Writing of the paper: H.F.L. and W.C.C.

BMD

bone mineral density

BP

biological process

eQTL

tissue expression quantitative trait loci

FDR

false discovery rate

FN

femoral neck

GEO

Gene Expression Omnibus

GO

Gene Ontology

ITPKC

Inositol-trisphosphate 3-kinase C

LD

linkage disequilibrium

LS

lumbar spine

SNP

single-nucleotide polymorphism

SOCE

store-operated calcium entry

STIM1

stromal interacting molecule 1