Integrative Analysis Confirmed the Association between Osteoprotegerin and Osteoporosis

Hui Tang, Xiaowei Zhu , Longfei Wu Xingbo MoFeiyan Deng Shufeng Lei＊

1Center for Genetic Epidemiology and Genomics, School of Public Health, Jiangsu Key Laboratory of Preventive and Translational Medicine for Geriatric Diseases,Medical College of Soochow University, Suzhou, Jiangsu 215123, China

2Zhangjiagang Center for Disease Control and Prevention,Zhangjiagang, Jiangsu 215600, China

Key words: osteoprotegerin; osteoporosis; single-nucleotide polymorphism;integrative analysis Objective This study aimed to verify the association between osteoprotegerin gene (OPG) and its variants with osteoporosis (OP) by performing integrative analysis.Methods We used the KGG software to perform gene-based association analysis, which integrated all publicly available single-nucleotide polymorphism (SNP)-based P values and obtained an overall P value for the OPG. The significant SNPs were screened for expression quantitative trait loci (eQTLs). Meta-analysis was used to combine the associations between the variants of OPG and bone mineral density (BMD) reported in the literatures. Then we performed dual-luciferase reporter gene systems for the functional verification of the variants of OPG in vitro.Results In the gene-based association analysis, the over all P value of OPG was 6.24×10-13 for BMD at femoral neck (FN)and 7.37×10-17 for BMD at lumbar spine (LS), indicating the importance of OPG for OP. The publicly available eQTL database identified 5 eQTLs which exert cis-regulation effects on OPG at FN and LS. Literature searching found that rs2073617 (known as T950C) was the hot spot SNP. There were 13 relevant studies on rs2073617 besides the GEFOS-2 study identified from the PubMed. Significant differences among TT, TC and CC genotypes at FN (P= 0.047) and LS (P=0.025) were shown by meta-analysis, demonstrating the associations between T950C polymorphism and BMD. Luciferase gene expression was significantly higher at the presence of allele C than allele T in the 293T cells (t=-9.47, P＜0.01).Conclusion The integrative analysis further confirmed the importance of OPG in OP and the correlation of T950C polymorphism with BMD of OP. The strategy can be used as a reference for functional interpretation of other disease-related genes.

O STEOPOROSIS (OP) is a systemic skeletal disease characterized by low bone mass and microarchitectural deterioration of bone tissue， with a consequent increase of bone fragility and susceptibility to fracture. OP is a global health threat in the elderly. Some studies estimated that the number of patients with osteoporotic hip fracture worldwide was more than 200 million，and 20% of hip fracture patients were complicated by pulmonary embolism and pneumonia， or died within 1 year.1

OP is under strong genetic determination， with the heritability of 60%-85%.2，3Previous genome-wide association studies (GWASs) have identified many genetic loci that are associated with OP.4-6However， the associations are largely inconsistent among different studies. Moreover， beyond the established statistical associations between genetic markers and OP at DNA level， functional mechanisms underlying the associations were rarely explored. Further functional experiments and more in-depth bioinformatics analyses are in urgent need to verify the importance of the identified genes and/or variants on OP.

Previous GWASs have consistently identified single nucleotide polymorphisms (SNPs) in osteoprotegerin gene (OPG) that are associated with OP.4-6OPG is located at 8q24. The product osteoprotegerin (OPG)，also known as TNFRSF11B (tumor necrosis factor receptor superfamily， member11b)， inhibits bone resorption mainly by OPG/RANK/RANKL axis where OPG plays a role as a soluble decoy for binding RANKL and weakens the inhibitory effect of RANKL on RANK activation. OPG is an important extracellular regulator for osteoclast development， inhibiting mature osteoclast formation and the survival of pre-existing osteoclasts.7，8

In this study， we integrate the publicly available results in the literature and datasets as well as our inhouse molecular experiments. The integrative analysis started with gene-based association analysis to evaluate the entire effect of OPG on OP; then we checked the results of expression quantitative trait locus (eQTL)analysis to test the regulation effects of the studied variants on mRNA expressions; subsequently， meta-analyses focused on SNP was performed to combine the association of existing results; finally， we performed luciferase reporter assays to validate the function of the focused SNP in vitro. The purpose of this study is to verify the association between OPG and OP and to identify certain SNPs within OPG that are closely related to OP.

MATERIALS AND METHODS

Osteoprotegerin gene-based association analysis

The data for the gene-based association analysis were acquired from a meta-analysis performed by the GEnetic Factors for OSteoporosis consortium (GEFOS-2)，9which investigated GWAS for bone mineral density(BMD) at femoral neck (FN， n=32，961) and lumbar spine (LS， n=31，800) based on 17 studies. The populations in these studies across North America， Europe，East Asia and Australia， with a variety of epidemiological designs and patient characteristics.9The BMD association P values of about 2.5 million SNPs at FN and LS were publicly available.10Data files for association analysis， including the rs number， chromosome， position and SNP-based association P value of each SNP for KGG were prepared using the R program. All the SNPs on OPG with BMD association P＜0.05 in the dataset were used for gene-based association analyses using the GATES (Gene-based Association Test using Extended Simes procedure) method11with KGG software.12We extended Simes test to combine P values based on all SNPs within one gene and obtained an overall P value for the gene.

Expression quantitative trait locus confirmation

Meta-analysis of SNP within OPG in relation to BMD

Meta-analysis study was performed in compliance with the STrengthening the REporting of Genetic Association studies (STREGA) statement， an extension of the Strengthening the Reporting of Observational Studies in Epidemiology (STROBE) statement.14Procedure of retrieving studies， extracting data and assessing quality of included studies were implemented by two authors independently and the discrepancies were discussed by all authors till reach consensus.

Retrieval of studies

We searched the MEDLINE database via PubMed in August 2017 using terms of “bone mineral density”， (“osteoprotegerin” OR “OPG” OR “TNFRSF11B”)， “polymorphism” and “humans”. The study subject was restricted to human. No language restriction was set. The potential eligible studies were screened by reading titles and abstracts. The full texts were available to eliminate the non-conformity studies. In addition， we used “related citations” function in PubMed and screened reference lists of all eligible studies to identify potentially eligible studies.

Inclusion criteria

Eligible studies that meet the following criteria were qualified for inclusion： (1) the mean and standard deviation (SD) values of BMD among the different genotypes of focused SNP were available; (2) the BMD was measured by dual-energy X-ray absorptiometry(DEXA); (3) the participants didn’t suffer from certain diseases or conditions that potentially affect bone development， such as chronic disorders involving vital organs， severe skeletal diseases， or malignant tumors;(4) the study must furnish enough data for further quantitative analyses， such as the number of subjects and BMD measurement sites.

Data extraction Two authors independently extracted the data using a standard form including the following items： (1) first author’s name， publication year and the geographic region where the study performed; (2) sample size， age，gender， menopausal status; (3) instrument of BMD measurement; (4) baseline of mean BMD and their corresponding SDs at different anatomic sites in each genotype group of focused SNP. The combined BMD of homozygote or heterozygote was also extracted， or calculated by combining the mean BMDs weighted by its proportion if the study did not contain the BMD values of three genotypes.

甲状腺乳头状癌为发病率较高的内分泌系统恶性肿瘤，分化好，恶性程度较低。其发病与性别、年龄、碘的摄入等多因素有显著的相关性[13-15]。随着对其研究的不断推进，我国甲状腺癌的诊治逐渐规范，甲状腺癌诊断方法由单纯使用超声发展为影像定位-细胞病理-分子靶标多层次早期诊断体系，诊断过程较为系统，具有科学性[16]。TPC恶性程度较低，主要治疗方案以外科手术为主。孟晓敏等对甲状腺乳头状癌原发灶外科治疗进行meta分析得出：TPC原发灶外科治疗术后生存率较为可观，且甲状腺部分切除术损伤喉返神经或甲状旁腺发生率较低[17]。

Quality assessment

Quality assessments were performed by using the agency for healthcare research and quality (AHRQ)assessment with 11 items for observational studies.15The full texts of the articles were carefully read and if answer for an item was “yes”， the item was scored 1; if“no” or “unclear”， then the item was scored 0. Included studies were divided into 3 categories： low quality(score 0-3); moderate quality (score 4-7); and high quality (score 8-11). All articles identified as low quality were excluded from the study.

Statistical analysis

The heterogeneity of the included studies was examined by using I-square. Fixed effect model was used to compute weighted mean difference (WMD) when there was no significant heterogeneity in the study(i.e. P＞0.1 and I2＜50%)， otherwise， the random effect model. The Z-test was performed to evaluate the significance of pooled WMD. Publication bias was appraised by Begg rank correlation test and Egger linear regression test. Sensitivity analysis was taken when there was significant publication bias. The meta-analysis was performed with STATA software (version 12.0，Stata Corporation， USA). Two-side P values were used and P＜0.05 was considered as statistically significant.

Luciferase reporter assays

Luciferase reporter assays were performed to validate whether the identified SNP had regulation effect on gene transcription. The pGL3 firefly luciferase reporter vectors (pGL3-enhancer vectors)， the pRL-TK renilla luciferase contrast vectors (pRL-TK vectors)， and dual-luciferase reporter assay system were purchased from Promega Corporation (Madison， WI， USA). The transfection reagent， lipofectamine 2000， was purchased from Invitrogen Corporation (Carlsbad， CA， USA).Plasmid purification kit was purchased from QIAGEN Corporation (Hilden， Germany). The transfected cells，293T cell line， were purchased from the cell bank of Chinese Academy of Sciences (Shanghai， China). 293T cells grew attaching to the wall and were fostered in high glucose media at 37˚C with a mixture gas of 5%CO2.

Target gene fragments， harboring wild and mutant alleles of the focused SNP were obtained by polymerase chain reaction (PCR). DNA samples were extracted from peripheral blood of a healthy human being to acquire wild fragment. Mutant fragment was obtained by point mutation. The PCR products containing the focused SNP were digested by KpnI and XhoI and linked into pGL3-enhancer vectors containing fluorescein gene to construct recombinant plasmids. Furthermore， sequencing of recombinant

plasmids confirmed the accuracy of target gene fragments. Both recombinant plasmids and pRL-TK vectors were amplified in Escherichia coli (E. coli) with Luria-Bertani (LB) media and purified with QIAGEN plasmid purification kit. Eventually， all plasmids for dual-luciferase reporter assays were sequenced and verified by PCR.

The assays were implemented using the dual-luciferase reporter gene systems following manufacturer’s instructions. The 293T cells were transfected with recombinant pGL3-Enhancer vectors along with pRLTK vectors by using lipofectamine 2000 respectively to verify the transcriptional activity of OPG. 293T cells were split and specific luciferase substrates were added. The luciferase assays can conveniently explore the efficiency of gene expression in mammalian cells with the firefly to renilla luciferase signal ratio， removing background signal value. Luciferase signals were measured using Enzyme-linked immunosorbent assay in Microplate Reader (BioTek， Winooski， USA). In the luciferase assays， we set three parallel holes for each group and take the means as results.

The difference of gene expressions between constructs with wild or mutant alleles was statistically compared using t-test by SAS software (version 12.0).We considered the test as two-sided and the differences among variants were significant when P＜0.05.

RESULTS

Gene-based association of OPG with OP

Association P values for 74 SNPs in OPG (from upstream 2-kb to downstream 2-kb) were analyzed，among which 62 SNPs for FN-BMD and 44 SNPs for LS-BMD with association P value ＜0.05 were identified. After combining the whole genetic information hinted by all SNPs in OPG， the gene-based P value was 6.24×10-13for FN-BMD and 7.37×10-17for LS-BMD，which confirmed the importance of OPG on OP.

EQTL confirmation and screen of focused SNP

We searched the corresponding eQTLs according to the identified 62 SNPs for FN-BMD and 44 SNPs for LSBMD， and found 5 SNPs exerting cis-regulation effects on OPG at the two skeletal sites with P＞0.01 (Table 1). The scores reported the rs2073617 (known as T950C， located in the promoter region)had the most different expression in different genotypes. Literature searching for the identified SNP found that rs2073617 has 13 relevant studies in the MEDLINE database apart from studies presented in GEFOS-2，indicating it was a noteworthy SNP.

Associations between T950C polymorphism and BMD

A total of 121 articles were retrieved according to the literature searching strategy. By reading the abstracts and the full texts， we excluded articles according our including criteria. Eventually， 7 articles were included for meta-analyses (Table 2).16-22These articles were published in the period from 2002 to 2009， and the subjects were from 6 countries. Data incorporated into our analyses were from 4，135 individuals.

Evaluation of BMD among subjects with different genotypes at both LS and FN showed that subjects with TT， TC and CC genotypes were significantly different in FN-BMD (P=0.047) and LS-BMD (P= 0.025).Fixed effect model was selected because of the insignificant heterogeneities (I2≤50% and P＞0.1) of the BMD in different genotypes in pairs. As shown in Figure 1 and Figure 2， BMDs at both LS and FN were significantly lower in subjects with TT genotypes than in subjects with TC (LS-WMD： -0.01; 95% CI：-0.01， -0.01; FN-WMD： -0.00; 95% CI： -0.00， -0.00)or CC (LS-WMD： -0.01; 95% CI： -0.02， -0.01; FN-WMD： -0.01; 95% CI： -0.01， -0.01) genotypes. Subjects with CC genotypes had significant higher BMD at LS (z=4.61， P＜0.001) and FN (z=12.78， P＜0.001)than subjects with the TC genotypes). Furthermore，subjects with TT homozygote had lower BMD and subjects with CC homozygote had higher BMD at both LS(z=13.00， P＜0.001) and FN (z=13.58， P＜0.001) than others.

Table 1. The eQTL P values from the publicly available dataset

Table 2. Basic characteristics of the studies included in the meta-analysis

Publication bias was not found in all models from the evidence of the results of Begg rank correlation test (P＞0.05) and Egger linear regression test (P＞0.05 and 95% CI including 0).

Function of T950C in vitro

Luciferase reporter assays validate the function of the identified T950C (Figure 3A) in vitro. Promoter fragments of 403bp harboring wild or mutant alleles were designed and sequenced (Figure 3B, 3C). After constructing the dual-luciferase reporter plasmids， we evaluated the effects of T950C polymorphisms on luciferase gene expression. As shown in Figure 3D， 293T cells with mutant allele T presented a quarter of luciferase signal than 293T cells with wild allele C (t=-9.47，P＜0.01)， indicating gene expression was significantly higher in the presence of allele C in the 293T cells.

DISCUSSION

Previous extensive GWAS studies have established the statistical associations between genetic variations and OP. However， it is a huge challenge to disclose functional mechanisms underlying the associations. Currently， such functional mechanisms are extremely complex and have not been fully explored.In order to cope with these complex mechanisms， it needs to integrate various aspects of evidence， especially those from publicly available data and results that archived in databases. Thus， integrative strategy on OP is in urgent need. It is effective by integrating multiple in-depth bioinformatics analyses with further functional experiments to confirm the importance of the identified genes and/or variants.

In this study， we combined the whole genetic information hinted by all SNPs in OPG. The gene-based P value confirmed the importance of OPG on OP. Genebased association analysis takes a gene as basic unit for association analysis. As this method can combine genetic information given by all the SNPs into a single gene， it can obtain more informative results and evaluate the entire effect of a gene on phenotype.23This method has been used as a novel complement method for SNP-based GWAS in identifying disease susceptibility genes.24-26

As genetic variants may alter binding affinity between variants and regulatory factors， and thereby alter transcription of mRNA and translation of target protein， i.e.， genetic variants may exert allelic specific influence on gene expression and play a crucial role in phenotypic variation or disease susceptibility. EQTL analysis is a powerful tool to evaluate the regulation effect of SNPs on mRNA expression. From the existing eQTL database， we found T950C had the largest variance in gene expression among different genotypes.

The meta-analysis of T950C and OP announced the differences in the associations of three genetypes with BMD. This result is consistent to the GEFOS-2 meta-analysis， where allele T was found to have negative correlation with LS-BMD and FN-BMD (P=6.19×10-15and P=2.50×10-9， respectively).

Finally， the relationship of the T950C polymorphism and OPG expression level were confirmed by the in-house functional experiments (luciferase reporter assays) in 293T cells. The result was identical with several functional experiments in HeLa， COS-7 and RAW 264.7 cells.27

OPG has been identified as an important candidate gene for OP. The frequently mentioned function of OPG is to act as decoy receptor for RANKL and neutralizes the effect of RANKL on osteoclastogenesis. Additionally， OPG has other important functions related to bone. For example， OPG also interacts with TNF-related apoptosis inducing ligand (TRAIL)and protect against the TRAIL-induced cell apoptosis.OPG-TRAIL binding blocks the inhibition of OPG-mediated osteoclastogenesis.28A lower OPG/RANKL ratio is associated with down-regulation of some important BMD associated genes (e.g.， RUNX2， COL1A1) in patients with osteoporosis.29Several other studies have discovered the significant correlation between BMD and a series of SNPs in the promoter region of OPG.Although these studies were not included in our meta-analysis because they did not meet our inclusion criteria， their results still have provided additional evidence to support association between OPG and OP. Moreover， these studies further investigated the in-depth functions of the detected variants. For example， multiple SNPs within OPG are associated with bone turnover markers， N-terminal propeptide of type I procollagen (PINP)， C-terminal cross-linked telopeptide of type I collagen (CTX-I)， etc.30Molecular haplotypes of OPG promoter region have differentially transcriptional action and show an apparent and exclusive interaction with a proximal SNP.31

Figure 2. Meta-analysis for the associations between T950C polymorphism and FN-BMD.FN-BMD： femoral neck-bone mineral density.

Extensively accumulating evidences have shown that the T950C polymorphism is functional variant with important regulation effect on OPG expression.32First， this polymorphism is located in the promoter region，16where DNA sequences can bind RNA polymerase or other transcription factors to control the mRNA transcription. For example， in the ChIP-Seq experiments，33the locus was bound with early B-cell factor 1 (EBF1)， which is a protein coding gene and recognizes variations of the palindromic sequence as transcriptional activator. Second， the relationship of the T950C polymorphism and OPG expression level was confirmed using eQTL analyses in a relevant study.34Third， the relationship of the T950C polymorphism and OPG expression level were confirmed by our in-house functional experiments and previous functional experiments.27All these experiments have consistently demonstrated that allele C has corresponding higher expression of OPG than allele T.

Figure 3. Diagram of the Dual-luciferase reporter assay.A. Relative position of T950C to the transcription initiation site of OPG. B and C， circle maps of dual-luciferase reporter gene vector containing mutant allele T or wild allele C. D. effects of T950C polymorphism on luciferase gene expression in the 293T cells. luc+， cDNA encoding the modified firefly luciferase; Ampr， gene conferring ampicillin resistance in E. coli; f1 ori， origin of replication derived from filamentous phage; ori， origin of plasmid replication in E. coli. Arrows within luc+ and the Ampr gene indicate the direction of transcription; the arrow in f1 ori indicates the direction of ssDNA strand synthesis.** t-test： t=-9.47， P＜0.01.

In summary， the results from our integrative analyses provide additional evidence supporting the importance of OPG and its variant (T950C) in OP by multiple in-depth bioinformatics analyses with further functional experiments (gene-based association，eQTL confirmation， meta-analysis for association， and function assay of SNP). It further improves our understanding in complex relationship between OPG and OP. The methodology in this study probably provides a reference for the analogous functional interpretation of other disease associated genes.

Acknowledgement

We would like to thank the Gilad/Pritchard group for their supply of the eQTL resources.

Conflict of interest statement

All authors declared no conflict of interests.