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

Osteoporosis is defined by low bone mineral density (BMD), which is mainly due to the imbalances in osteoclast and osteoblast activity. Previous studies indicated that early activation of osteoclasts relies on calcium entry through store-operated calcium (SOC) entry, and several genes, including STIM1, ORAI1, and ITPKC, are known as key regulators of SOC entry. However, the relationships between STIM1, ORAI1, ITPKC, and human BMD are still unclear. In order to investigate the plausible associations between these genes and BMD, we conducted a meta-analysis of genes expression and BMD using the publicly available GEO database. We further recruited 1044 subjects and tested associations between polymorphisms in these genes and BMD. Clinical information (including age, sex, and BMI) was collected and used for the analysis. Our results indicated that ITPKC gene expression was significantly associated with BMD. Furthermore, we found that one ITPKC SNP (rs2607420) was significantly associated with lumbar spine BMD. Through bioinformatics analysis, rs2607420 was found to be very likely to participate in the regulation of ITPKC expression. Our findings suggest that ITPKC is a susceptibility gene for BMD, and rs2607420 may play an important role in the regulation of this gene.


Introduction
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][13][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 3 ) [16,17]. Once IP 3 is phosphorylated to IP 4 , IP 3 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 3 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.

Dataset collection
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.

Meta-analysis of gene expression
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 2 . 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.

Patients and methods
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

Functional annotation of SNPs
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.

Statistical analysis
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.

Studies included in the meta-analysis
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.   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 2 = 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.

Demographic and clinical characteristics of subjects
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.

Haplotype associations for ITPKC and BMD at the LS and FN
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). The p value was adjusted for age, sex, and the body-mass index. p and q values of < 0.05 are shown in bold. q values of < 0.05 were considered statistically significant after correction for multiple testing.

Discussion
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 2 > 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.

Conclusion
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.