Glutamine Repeat Variants in Human RUNX2 Associated with Decreased Femoral Neck BMD, Broadband Ultrasound Attenuation and Target Gene Transactivation

RUNX2 is an essential transcription factor required for skeletal development and cartilage formation. Haploinsufficiency of RUNX2 leads to cleidocranial displaysia (CCD) a skeletal disorder characterised by gross dysgenesis of bones particularly those derived from intramembranous bone formation. A notable feature of the RUNX2 protein is the polyglutamine and polyalanine (23Q/17A) domain coded by a repeat sequence. Since none of the known mutations causing CCD characterised to date map in the glutamine repeat region, we hypothesised that Q-repeat mutations may be related to a more subtle bone phenotype. We screened subjects derived from four normal populations for Q-repeat variants. A total of 22 subjects were identified who were heterozygous for a wild type allele and a Q-repeat variant allele: (15Q, 16Q, 18Q and 30Q). Although not every subject had data for all measures, Q-repeat variants had a significant deficit in BMD with an average decrease of 0.7SD measured over 12 BMD-related parameters (p = 0.005). Femoral neck BMD was measured in all subjects (−0.6SD, p = 0.0007). The transactivation function of RUNX2 was determined for 16Q and 30Q alleles using a reporter gene assay. 16Q and 30Q alleles displayed significantly lower transactivation function compared to wild type (23Q). Our analysis has identified novel Q-repeat mutations that occur at a collective frequency of about 0.4%. These mutations significantly alter BMD and display impaired transactivation function, introducing a new class of functionally relevant RUNX2 mutants.


Introduction
Osteoporosis is the major metabolic bone disease among developed nations. The disease is characterized by low bone density and reduced bone quality through the deterioration of bone micro architecture. As a consequence, sufferers have impaired skeletal strength and an increased susceptibility to osteoporotic fractures. Bone mineral density (BMD) is a complex trait and is controlled by environmental and genetic factors [1]. BMD is a primary predictor of osteoporotic fractures; it is however a continuous trait related to age and weight. Osteoporosis as a category is a truncation of this continuous normally distributed trait, defined by a BMD value less than 2.5 SD from the population mean of young adults [2]. The proportion of fractures attributable to osteoporosis ranges from 10%-44% [3].
RUNX2 is a key regulator of skeletogenesis, bone and cartilage formation [4][5][6][7][8] and has been genetically associated with BMD [9][10][11][12][13][14][15]. RUNX2 transactivates genes such as osteocalcin, type 1 collagen and osteopontin and is initially expressed in mesenchymal condensations where it plays an essential role in osteoblast differentiation and in chondrocyte maturation [16]. RUNX2 deficient mice have un-mineralized skeletons with few immature osteoblasts and an absence of vascular and mesenchymal cell invasion in the cartilaginous templates. Heterozygous mutations in coding and promoter sequences of RUNX2 cause the dominantly inherited skeletal syndrome cleidocranial dysplasia (CCD). The disorder is characterised by persistently open or delayed closure of sutures, hypoplasia/aplasia of clavicles, Wormian bones, supernumerary teeth, short stature and other skeletal abnormalities [17,18]. A spectrum of severity of CCD symptoms is correlated with the extent of transcriptional activity remaining in the cognate mutant RUNX2 with severest symptoms associated with complete loss of RUNX2 function [19] and milder symptoms correlated with some residual RUNX2 transactivation function [20]. RUNX2 has consecutive polyglutamine and polyalanine repeats (Q/A repeat) in the protein sequence. Such repeat regions have the capacity to mutate via strand slippage during DNA replication. Glutamine repeat sequence expansion has been the cause of some diseases that show genetic anticipation, where severity increases in subsequent generations as repeat length increases due to errors in replication [21]. Wild type human RUNX2 contains a 23Q/17A repeat; 23 consecutive glutamine followed by 17 alanine residues. An insertion of the polyalanine tract (23Q/27A) was observed in one patient with CCD [7] although no further polyalanine-related CCD patients have been reported and no evidence currently associates the Q-repeat region with CCD. We previously identified Q-variants (15Q, 16Q, 24Q and 30Q) in an Australian fracture cohort [9] and two 16Q variants as well as a single alanine expansion variant (23Q/23A) in a randomly selected population from Aberdeen [10]. Pineda et al. [22] observed three Q-variants (16Q, 18Q and 30Q) in a Spanish population. Evidence exists that the RUNX2 Q-repeat is a site of functional variation; in carnivores the length of the Q repeat is significantly associated with mid-face length and nose curvature [23], indicating an effect on bone growth rates. Sears et al. [24] confirmed the functional nature of carnivore related Q/A repeat alterations, showing correlation of transcriptional variation with facial length in carnivores. These lines of evidence suggest that RUNX2 repeat variation may be functionally different in the human population. We hypothesised that glutamine repeat variants would exist in normal populations and may influence adult BMD and/or risk of fracture.

Study Subjects
The study subjects were obtained from four different epidemiological studies of bone density and are summarised in Table 1. All subjects were female. Voting in elections in Australia is compulsory and each region maintains an electoral roll. The West Australian study [25] consisted of females between the ages of 70 and 85 years who were randomly ascertained from the electoral role and approached by letter. The Geelong Osteoporosis study comprised a cohort recruited at random from the electoral role (GOS random) [26] and a specific fracture study (GOS fracture), where all fracture cases over the age of 35 years presenting at the only two radiology practices in the region were invited to join, as described [27]. Genotype was obtained on 834 subjects with no history of fracture from the random sample of GOS and 578 of the GOS fracture cases. The Sydney study comprised monozygotic and dizygotic twins who volunteered for studies of general medical issues through the Australian twin registry, as described [28]. Self reported history of fracture was available among a series of questionnaires completed by volunteers. Data were obtained on 980 participants. In the Tasmanian older adult cohort (TASOAC) study, 388 DNA of females were available of which 385 yielded genotype data. In this study, subjects ranged from 50 to 80 and were recruited by random ascertainment from the electoral roll [29]. For the Aberdeen study, subjects were postmenopausal women aged between 45 and 55 years who were approached at random using Community Health Index records of Aberdeen, Scotland as previously described [30]. Data were obtained on 991 study participants. In addition, 101 elderly female clinic patients with established osteoporosis and at least one osteoporotic fracture were available from the Geelong endocrine clinic. A similar cohort of 200 clinical patients with established osteoporosis defined by recurrent vertebral fracture was available from a New Zealand study, as described [31,32]. BMD data were not available for clinic patients. Individuals in all studies were female Caucasians and no subjects were excluded. Appropriate written consent was obtained from subjects under procedures in accordance with the Declaration of Helsinki approved by the relevant human research ethics committees (HREC) as previously described [25][26][27][28][29][30][31][32]. The committees are as follows: Sir Charles Gairdner Hospital HREC, Barwon Health HREC, Northern Sydney Local Health District HREC, Southern Tasmanian Health and Medical HREC and University of Otago HREC.

Bone measures
Bone density was measured by dual energy X-ray absorptiometry (DEXA). In the Aberdeen study a Norland XR26 or XR36 was used (Norland Corporation, Fort Atkinson, WI, USA). A Lunar DPX-L machine was used for the GOS. The Western Australian study used a Hologic 4500A fan beam densitometer to measure total BMD of the hip region. A Lunar Achilles ultrasound machine was used to measure the left calcaneus bone. Measures of speed of sound (SOS), broadband ultrasound attenuation (BUA) and bone stiffness were available. In the Sydney study, a Hologic 4500A fan beam densitometer was used to measure BMD at the lumbar spine, total hip, femoral neck and whole body. Western Australian, TASOAC and Sydney studies used the same type of densitometer and were calibrated with the same standard phantoms.

Detection of Q-variants
PCR was used to amplify RUNX2 exon 1 fragments harboring the Q/A repeat using the primers forward 59-CCGGCAAAAT-GAGCGACG-39 and reverse 59-GGGCGGTGTAGCCTCT-TACCTT-39. The PCRs were carried out in 20 ml reactions containing 561027 M primers, PCR reaction buffer, 125 mM dNTPs, 80 ng genomic DNA, 0.5 units Taq polymerase in 20 ml as specified by the supplier (Promega, Sydney, Australia). The study populations were genotyped separately by differing means: testing 100 random DNA by all methods gave the same genotypes. For the GOS and TASOAC, the resulting 336 bp fragments were resolved via nondenaturing 10% polyacrylamide gel electrophoresis (PAGE). Heteroduplex analysis was used to determine the presence of Q-repeat variants. For the Aberdeen populations the PCR fragment was digested with MspA1I (New England Biolabs) and resolved on 3% agarose. For the Western Australian samples, all PCR products were initially analysed via dHPLC (Varian Prostar 430, Varian Industries, Sydney, Australia) and subsequently by heteroduplex PAGE. Genotyping was done twice for all samples and all variant genotypes confirmed by PAGE to reveal heteroduplexes. All PCR amplified DNA fragments that displayed unique mobility patterns on PAGE were sequenced using BigDye

Plasmids, transfection and cell culture
Expression plasmids containing mutant Q/A domains were constructed by cloning the XbaI-XbaI fragment from full length RUNX2 expression plasmid pEF-aA [33] into pUC18 to create pUC18RUNX2 which served as a template to create glutamine variants via PCR cloning. PCR was used to amplify partial RUNX2 promoter and exon 1 fragments from DNA samples harboring Q/A mutations using the oligonucleotide primers 59-TTCACCACCGGACTCCAACT-39 for the 59 side and 59-CATCTGGTACCTCTCCGAGGGCTACCACCTTGAAGG-CCACGGGCAGGGTC-39 for the 39 side. The reverse primer contained an EcoNI tag facilitating the cloning of the PCR amplified product into the BglII-EcoNI site of pUC18RUNX2. Mutant Q/A regions were confirmed by DNA sequencing and the XbaI-XbaI fragments were cloned into the XbaI-XbaI site of the mammalian expression vector pEF-BOS [34]. Restriction digest analysis and DNA sequencing confirmed the orientations of RUNX2 inserts. pGL3-Basic served as a template to create the reporter constructs. BglII and HindIII sites were introduced by PCR amplification in the human osteocalcin gene using oligonucleotide primers 59-CAGGAGATCTCTGACCGTCGAGCTG-39 for the 59 side and 59-GGGCAAGCTTGGTGTCTCGGGTGG-C-39 for the 39 side. The resulting BglII-HindIII fragment was cloned into the corresponding sites of pGL3-Basic to create pOSLUC, which has 590 basepairs of the human osteocalcin promoter, including 70 base pairs of the untranslated message. Oligonucleotides containing a consensus mouse osteocalcin RUNX2 response element (RRE: 59-GATCCGCTGCAATCACCAACCACAGCA, with RUNX2 consensus binding site underlined) were cloned into the BglII site of pOSLUC to create pRRE, which had three copies of the RRE inserted as direct repeats. Transfections were performed using FuGENE (Roche) as the manufacturer's protocols. Antibody immuno-staining was performed as described in Yoshida et al. [35]. 1,25 dihydroxyvitamin D 3 (1,25(OH) 2 D 3 ) was kept in the dark under argon in isopropanol and added with appropriate vehicle controls after dilution in ethanol. No more than 0.1 ml ethanol per ml medium was present in cell culture. Western blots of transfected cells were done by standard methods using anti-RUNX2 antibody (D130-3, MBL International) according to the manufacturer's instructions.

Luciferase Assay
NIH3T3 cells and HEK293 cells were maintained in DMEM (Gibco) supplemented with 10% FBS (v/v) (Gibco), 1% Penicillin-Streptomycin (v/v) (Gibco) in a 5% CO 2 humidified atmosphere at 37uC. For the luciferase assay NIH3T3 cells were seeded into 6well plates at a density of 1610 5 cells/well 24 hours prior to transfection. HEK293 cells were seeded into 12-well plates at a density of 1610 5 cells/well 24 hours prior to transfection. NIH3T3 and HEK293 cells were transfected using FuGENE6 transfection reagent according to the manufacturer's instructions (Promega). Cells were harvested 48 hours post-transfection. Luciferase activities were determined using the dual luciferase assay system described by Dyer et al. 2000 [36]. pRL-CMV (Promega) was used as an internal control to normalize results.
Gel shift assays, in vitro protein translation and GST pulldown assays Gel shift assays were done by standard methods as previously described, using 5% PAGE and in vitro translated proteins [37]. DNA probes representing RUNX2 binding sites were end labelled using a fill in reaction of 59 overhangs using DNA polymerase 1 Klenow fragment and fragments purified by PAGE. In vitro translated human RUNX2 variant proteins 23Q wild type, 16Qvariant, and 30Q-variant were generated from DNA fragments representing 23Q, 16Q, or 30Q variants recloned into a T7 promoter containing hRUNX2 vector by using T7 RNA polymerase (TNTH Promega, USA) and coupled transcriptiontranslation reactions containing [ 35 S]-methionine. Protein products were visualised and quantified by [ 35 S]-methionine incorporation. Bacterial over-expression of glutathione-S-transferase coupled human VDR (GSThVDR) fusion protein was performed by induction with isopropyl-b-D-thio-galactopyranoside (IPTG, 1.25 mM) for 3 h at 30uC in the JM109 E.coli strain as described previously [37]. GST  GSThVDR-Sepharose was released from the Sepharose, electrophoresed through a 12% SDS-PAGE, and detected by autoradiography using either phosphor screens (BioRad) or X-ray film. Densitometric analysis on resolved bands was performed using the ChemicDoc system (BioRad).

Statistical methods
Within studies, Q-variants were analysed using analysis of variance (ANOVA) and Student's T-test. Genotype status (Qvariant carrier or not) was coded as one and zero, respectively, and used as a variable in ANCOVA analysis to test the effect of covariates such as age and weight. Age, age-weight and ageweight 21 adjusted values of bone density parameters were produced using linear regression. Stiffness was log e transformed to comply with equality of variance assumptions in ANOVA. These adjustments had no material influence on the conclusions and Z scores of age-weight 21 regression are presented. Incident fracture during five years of observation was available for the Western Australian study. Incident fracture was categorised as zero and one for absence and presence of fracture, respectively, for analysis by logistic regression and the number of such fracture events also examined in genotype groups using contingency tables. Results of logistic regression are presented as p values and odds ratios (OR) with 95% confidence intervals (CI). Genotypes were examined pooled (all variants) and separately as either deletions or insertions. The only measure in common across all studies was femoral neck BMD. In order to pool all studies, each individual was expressed within a study as a Z-score of residual femoral neck BMD from the age-weight 21 regression specific for that study. Under the null hypothesis of no effect at this locus on BMD, the Z scores should be distributed randomly around a mean of zero. Population simulations were constructed in excel, using the Monte Carlo simulation plugin, PopTools [38], and the correlated variable tool to simulate bone density parameters as a series of correlated normal distributions, based on the correlations observed in the populations studied. Results of transfections were analysed using ANOVA with Fisher's least significant difference test for pair wise comparisons. Counts of cells were analysed using contingency tables and Chi-square. Analysis of covariance was used to analyse categorical and continuous covariates.

Q-variants identified
From the 1078 individuals genotyped in the Western Australian study, there were five 16Q alleles, one 18Q allele and two 30Q alleles. In the GOS random population sample three Q-variants were identified in 822 subjects: these were one 15Q and two 16Q variants. In the GOS fracture study, where recruitment was based on any fracture at any age, four Q-variants were identified in 598 subjects: these were two 16Q and two 30Q variants. In a prior publication [9], three of these variants were reported (two 16Q and one 30Q): in this study more samples from the GOS fracture cohort were genotyped and an additional variant identified. In the Sydney population of 980, five individuals with Q-variants were identified: carrying 16Q, 18Q and 30Q. The 16Q and 18Q variants were found in two pairs of monozygotic twins. The 30Q variant was in one individual of a dizygotic twin pair. For analysis, the average bone density values of the monozygotic twin pairs were taken and the pair considered as a single genetically unique individual. In the Tasmanian study (TASOAC) no Q-variants were found in 385 genotyped. In the Aberdeen study, two 16Q variants were identified previously from 991 subjects [10]; in this study their bone parameters are reported. In order to estimate the frequency of the Q-variants without bias, only those populations where the recruitment strategy involved volunteers were considered (Western Australia, Geelong, Tasmania and Aberdeen). In this case, 13 Q-variants were observed within 3276 subjects, yielding a population frequency estimate of 0.004 (95%CI, 0.002 to 0.007). Using the binomial theorem and the observed frequency of Q-variants, it is estimated that 95% of repeat studies of a similar size (approximately 3000 subjects) will detect between 6 and 19 such Q-variant carriers, if allele frequency is similar in other populations. In 103 female post-menopausal clinic patients with established osteoporosis with fracture, 16Q and 30Q variants were found. In a similar population of vertebral fracture cases from New Zealand [31,32], two 30Q variants were found in 200 patients genotyped. Therefore, in clinic patients with osteoporosis, four Qvariants were found in 303 patients genotyped. Apart from those Q-variants described above, three examples of 24Q alleles were detected overall. These 24Q subjects were not considered in this analysis. All Q-repeat variants were heterozygous carriers with the other chromosome containing a 23Q/17A wild type allele.

Bone density in Q-variants
Characteristics of Q-variant carriers with respect to femoral neck (FN) Z scores are presented in Table 2. Figure 1 shows the mean Z score of Q-variants for each measured parameter in each study cohort. The Q-variant carriers in Western Australian were compared to non carriers (N = 1021) for differences in the ultrasound measures of bone density: broadband ultrasound attenuation (BUA) and speed of sound (SOS). Q-variants had significantly lower BUA measures (n = 8, BUA, p = 0.025) at the calcaneus, with a difference of 20.80 standard deviations (SD, or Z score) and significantly lower SOS (p = 0.05) and log stiffness (p = 0.013). Weight, height and age were not significantly different in Q-variants. For all BMD values, age and age-weight 21 adjustment was done using regression and adjusted values produced. Individual values were then expressed as Z scores around the relevant population mean. While the particulars of reported measures change depending on such adjustment, the general outcomes reported here are not dependent on this mathematical adjustment. Measures of bone density using DEXA were similarly lower in Q-variants (Fig. 1A) although only nominally significant with the total hip BMD measure (p = 0.028). At other sites there was a trend for Q-variants to have lower BMD compared with all others in the study: femoral neck (p = 0.22), trochanter (p = 0.065), and intertrochanteric region (p = 0.06). Of the seven Q-variants identified from GOS there was a trend for negative Z scores (BMD below the population mean value, Fig. 1A) with 6 of 7 sites having negative average BMD. The Q-variant subjects from the Aberdeen study had negative Z-scores for both femoral neck and lumbar spine BMD (averaging 21.39SD and 21.36SD, p = 0.01 and p = 0.11 respectively, Fig. 1A). In the Sydney study, the data from each monozygotic twin pair were averaged and thus three were available with BMD data. The Q-variants had negative BMD Z-scores for the femoral neck with an average effect of 20.62SD. There was a trend for negative Z scores for all bone measures (Table 3).

Combining BMD data using Monte Carlo simulation
Bone density measures at different anatomical sites are correlated: additional measures in individuals therefore contribute useful information to explore the bone phenotype of Q-variants. All these data from skeletal sites can be used collectively to explore the genotype-phenotype relationship. Empirical p values based on Monte Carlo modelling make use of the additional data available in bone studies, where the scanning methods provide BMD measures of multiple sites. Considering only Z scores derived from within each study eliminates the effect of using different DEXA devices. A simulation was made based on the likelihood of sampling such individuals with measurements of correlated bone related variables. Correlated variables were simulated based on the measured correlations of bone parameters in the various populations, and a simulation of sampling was done in order to determine how many such samples met or exceeded the observed data. This approach was first taken with the Western Australian study, considering just deletion variants. There were six deletion variants, identified in a population sample of 1078 subjects, each with seven measured parameters related to bone. In this case, of 100,000 simulations of sampling six individuals measured for seven correlated bone parameters, only 252 such samples had an overall group mean equal to or more negative than the observed mean of 20.74 Z scores for all variables. This yields an empirical multivariate p value of 0.0025. The GOS study came from a similar population via similar recruitment with seven bone related parameters measured; although the bone parameters were measured using different equipment. The three deletions observed in GOS study were combined with those from Western Australia and a simulation done based on the idea of sampling from the combined population: in this case the observed mean was 20.63 Z scores and 235 samples from 100,000 trials met or exceeded the observed mean, giving an empirical p value of 0.0024. When 30Q variants were combined with all other variants from the two Australian studies the Monte Carlo based multivariate p value was: 0.0007. Of the 75 bone related measures from the Q variants identified in Aberdeen, Western Australia and Geelong studies, 59 measures had negative Z score values. Combining data from the three studies for Monte Carlo simulation resulted in 9 samples in 100,000 trials; giving an empirical probability of p = 0.0009. These simulations suggest that the additional measures provide useful data supporting the hypothesis that Q-variant carriers have lower bone density and that the average deficit is of the order of 0.7 Z scores of age, weight 21 adjusted BMD.
Using BMD data from the Q-variants derived from the GOS fracture cohort is not necessary, since a significant p value was already obtained using data from studies with no apparent recruitment bias. However, if the GOS fracture and Sydney cohort data were used in the Monte Carlo simulations the result was still significant. Addition of all subjects for whom BMD data were available to the Monte Carlo simulation supported the hypothesis that Q-variant carriers have significantly lower BMD (3 samples in one million trials with mean values equal to or more extreme than the observed means).
Femoral neck BMD was the only measure in common between the GOS random sample, GOS fracture, Aberdeen, Sydney and Western Australian studies. In order to combine the effects of these studies, age-weight 21 adjusted Z-scores for each study were taken giving a mean Z score of 20.60. Simply combining all populations, as if there were one global population from which the samples were drawn, gave a p value of 0.0007. The conclusions were not altered appreciably by taking the average of those who were twins.

Relationship to fracture
In the GOS study, the random sample that was genotyped excluded any with history of fracture, whereas in the Western Australian study, although recruitment was essentially random, fracture was not a reason for exclusion. Of the eight Q-variants found in the Western Australian study, four had prior osteoporosis related fractures (hip, spine and arm fractures). Data on incident fracture during a five year observation period were available for those who had completed the follow-up. Incident fracture during the study period was coded as 0 or 1 for absence and presence of fracture and the type of fracture was ignored. Four of 8 Q-variant carriers sustained incident fracture whereas the fracture rate was 17.8% in 1036 subjects with the normal 23Q RUNX2 allele. Qvariant carriers were significantly more likely (Fisher's exact test, UD, ultradistal radius BMD. Error bars are standard error. Z-scores presented are the within-study age-weight 21 adjusted Z score derived from within that cohort. B and C. Q-variant expression vectors have diminished capacity to transactivate target promoters in NIH3T3 mouse fibroblasts (B) and human embryonic kidney cell (HEK293) (C). Cells were transfected as described in the methods. In B the target was pRRE which is the human osteocalcin promoter with three additional synthetic mouse RUNX2 target sites added upstream, and in C the target promoter construct was 590 bp of the authentic human osteocalcin promoter driving luciferase (pOSLUC). D. Transfected cells show appropriate sizes of expressed product and exhibit equivalent amount of protein. NIH3T3 cells were transfected with empty vector (Con) or expression constructs as indicated. Cellular extracts were Western blotted after electrophoresis on SDS-PAGE. Numbers indicate molecular weight markers in kilodaltons. In B and C transactivation is expressed as percentage relative to the 23Q variant in order to normalise scales of NIH3T3 and HEK293 cells. doi:10.1371/journal.pone.0042617.g001 p = 0.036) to be in the fracture category (OR 4.7 with 95% CI 1.2 to 19.2). Overall, 6 of 8 Q-variants from the Western Australian study had prior or incident fractures. If the Geelong fracture study (n = 598) is combined with Western Australian study subjects who had ever fractured (n = 426), 1024 Australian subjects are available who had sustained fracture. Combining Western Australian [25] and GOS [39] subjects who had not reported fracture gave 1478 Australian non-fracture subjects. Despite the fact that no age matching was done, there was a significant increase (p = 0.036) in Q-variants within the fracture category in this comparison. Overall, of 22 Q-repeat variants identified, 12 had sustained some form of bone fracture, including the clinically important sites of femur, hip, spine and arm. Within the 30Q variants, three of five had fractured.

Reduced capacity to transactivate RUNX2 target genes
Expression vectors were constructed to express 23Q, 16Q and 30Q variants of RUNX2. Significantly lower transactivation of target gene promoters occurred after transfection for 16Q and 30Q compared to the wild type 23Q was observed using a semisynthetic RUNX2 reporter construct (pRRE) in mouse fibroblastic NIH3T3 cells (p = 0.002 and 0.016, for 16Q and 30Q, respectively, Fig. 1B). Similar data was obtained using the human osteocalcin promoter driving luciferase (pOSLUC) in HEK293 cells (Fig. 1C). Western blot of total protein from transfected NIH3T3 cells demonstrated the expected differences in apparent molecular weight in PAGE while there was no evidence a difference in relative abundance of Q-variants compared to 23Q wild type. Control cells transfected with empty vector showed no immunoreactivity (Fig. 1D).

Comparison with CCD-related RUNX2 mutants
CCD is a disease where severe inactivating mutations of RUNX2 exist. Some CCD-related RUNX2 mutants have no transactivation capacity while others, related to milder phenotypes, have a residual activity. CBFB is a heterodimer binding partner of RUNX2, known to increase DNA binding. The ability of RUNX2 variants to co-localise with endogenous CBFB in either nuclear, and/or cytoplasmic compartments was assessed by confocal immunomicroscopy with simultaneous staining with RUNX2 and CBFB with different fluorophores (Fig. 2 A-C). Cells were counted with exclusively nuclear staining, nuclear plus cytoplasmic and cytoplasmic only staining. The data were expressed as percentages of cells counted with a particular type of staining. For RUNX2, more cells with cytoplasmic staining were observed in transfected cells with 16Q and 30Q variants than expected compared with wild type 23Q (p = 0.007), with the 30Q having the stronger effect (p = 0.01). The majority (90%) of cells transfected with wild type 23Q RUNX2 also showed nuclear localization of CBFB: the other 10% represented cells with both nuclear and cytoplasmic staining. For transfected 16Q variant, 75% of positive cells showed nuclear RUNX2 while only 45% showed nuclear CBFB (p = 7.5610 25 ) with the difference appearing in the combined nuclear and cytoplasmic compartment (Fig. 2C). For transfected 30Q, the difference in CBFB binding was not as great: 64% showed nuclear RUNX2 while 54% of cells were also being positive for nuclear CBFB (p = 0.002) with a 10% increase in the nuclear and cytoplasmic compartment for CBFB.
The transactivation functions of the Q-variants were analysed using a third RUNX2 target gene reporter assay and activity compared to RUNX2 mutants known to cause CCD. Using the mouse osteocalcin reporter (p147) in NIH3T3 mouse fibroblastic cells: 16Q and 30Q variants had reduced capacity to transactivate  and this was a stronger effect than had been observed in the prior reporter-cell line combinations (Fig. 2D). The depressed transactivation ability of the 30Q variant was approximately the same as that observed in the Q266X CCD mutant, whereas the 16Q variant was similar to that seen in the A348fs mutant. Exogenous transfected CBFB enhanced wild type RUNX2 activity and rescued, to a large extent, the transactivation activity on the p147 target promoters when Q-variants were considered (Fig. 2E). Under the same conditions, the activity of CCD-related RUNX2 mutants was not rescued fully by co-transfection of CBFB partner (Fig. 2E). Endogenous CBFB nuclear translocation was significantly related to RUNX2 Q-variants when observed using confocal microscopy ( Fig. 2 A-C). Despite the fact that endogenous CBFB localization might be altered in RUNX2 variants, we reasoned that transfection of CBFB expression vector, may alter the phenotype of the Q-variants by augmenting endogenous CBFB levels. Expression vector for CBFB was transfected with Q-variants and activity on the p147 target in NIH3T3 cells examined.

Vitamin D receptor (VDR) and Q variants
In order to test the interaction of VDR with RUNX2 Qvariants in cells, the human osteocalcin promoter luciferase construct (pOSLUC) was transfected into NIH3T3 fibroblasts, with or without VDR expression vector and with or without 1,25(OH) 2 D 3 and RUNX2 expression vector (Fig. 3A). Differences in target promoter activity were once again observed, as described above, with 16Q and 30Q showing significantly lower activity compared to wildtype 23Q (Fig. 3A). As expected, transfection of VDR resulted in induction of pOSLUC. In the presence of 45 ng RUNX2 vector, 10 ng VDR vector and 1,25(OH) 2 D 3 (at 10 28 M) approximately 63 fold induction of the target vector occurred compared to empty vector vehicle treated cells (Fig. 3A). In the absence of 1,25(OH) 2 D 3 , transfected VDR had no significant effect on promoter activity. Treatment of cells with 1,25(OH) 2 D 3 , in the absence of transfected VDR, resulted in some induction of pOSLUC, consistent with activation of endogenous VDR. Induction by 1,25(OH) 2 D 3 was increased in the presence of transfected VDR (p = 1.7610 210 ). In the absence of 1,25(OH) 2 D 3 , both 16Q and 30Q variants had significantly lower target gene RUNX2 Q-variants compared with CCD-associated RUNX2 mutants, tested against the p147 mouse osteocalcin construct. In this assay Q-variants 16Q and 30Q have a deficit of transcriptional activity that overlaps that of CCD-associated variants. E. Co-transfection of CBFB expression vector with RUNX2 Q-variants completely overcomes the defect in transcriptional activity detected by the p147 promoter. In contrast, the transcriptional activity of CCD-associated RUNX2 mutants on the p147 construct is not rescued by CBFB co-transfection. The transcriptional activities of 16Q and 30Q alleles were compared to wild type 23Q RUNX2 and known CCD-related mutants using an in vitro target gene reporter assay consisting of the proximal mouse osteocalcin promoter (p147) driving luciferase target gene construct and transfected with expression vector for RUNX2 variants in mouse NIH3T3 fibroblast-like cells. In D and E, transactivation data are not percentages relative to 23Q wild type, but arbitrary units related to the ratio of firefly luciferase to Renilla luciferase. Key to symbols: Con, control cells with target gene promoter p147 but transfected with empty expression vector; 23Q, wild type; 16Q and 30Q, Q variants; R211Q, missense mutation at residue 211; Q266X, early termination at residue 266; A348fs, frames shift mutation at alanine residue 348. doi:10.1371/journal.pone.0042617.g002 promoter activity compared to 23Q wild type construct (23Q versus 16Q, p = 1.1610 25 ; 23Q versus 30Q, p = 4.0610 26 ), regardless of the presence or not of transfected VDR. In contrast, in the presence of 1,25(OH) 2 D 3 and transfected VDR the difference between 23Q and other Q-variant constructs measured by osteocalcin promoter activity was eliminated or reduced (23Q versus 16Q, p = 0.95; 23Q versus 30Q, p = 0.08).
VDR is a RUNX2 partner protein that can be activated by ligand, 1,25(OH) 2 D 3 . To examine the interaction of VDR and RUNX2 variants, in vitro translated proteins were made and tested in glutathione-S transferase (GST) pull down reactions for differences in protein-protein interaction. Unlabelled GST-VDR fusion construct was used to pull down 35 S labelled Q variant RUNX2 proteins (Fig. 3B). Under these conditions, GST-VDR demonstrated a strong interaction with RUNX2 (approximately 25% of label being recovered in pull down) but no significant difference in capacity to pull down RUNX2 Q-variants (p = 0.08), suggesting that protein-protein interactions of RUNX2 and VDR are unaffected by either variation in the glutamine repeat or the GST moiety added to VDR. There was no difference in the Figure 3. RUNX2 Q-variant interaction with the vitamin receptor (VDR). A. Transfected VDR, activated by 1,25(OH) 2 D 3 , can overcome the deficit in Q-variant transcriptional activity on a target human osteocalcin promoter. Although 16Q and 30Q variants show significantly reduced transactivation in NIH3T3 cells, as seen previously, transfected VDR alone does not alter the significant defect, however activation of VDR by the natural ligand 1,25(OH) 2 D 3 , results in essentially no difference in activation on the target gene promoter between the Q-variants. Promoter activity refers to human osteocalcin promoter construct driving fire fly luciferase (pOSLUC). Transactivation units are arbitrary based on the ratio of fire fly to Renilla luciferase. Empty vector controls are not shown but had activity of about one half unit on this graph. B. There was no detectable difference in the affinity of [ 35 S] RUNX2 Q-variants for binding labelled VDR protein in in vitro pull-down assays where VDR was used to pull down 35 S labelled RUNX2 protein via an affinity column. C. The ability of in vitro translated RUNX2 Q-variants to bind a RUNX2 DNA element from the mouse osteocalcin promoter (known as ORE) is enhanced by VDR in a ligand dependant manner. DNA binding was measured by phosphor analysis of labelled bands in gel shift assays after 5% PAGE. A difference in binding capacity existed on this element. Although addition of VDR or VDR and 1,25(OH) 2 D 3 , resulted in a increased total binding, in a ligand dependent manner, there was still a difference in Q-variants relative to 23Q. D. Images show gel shift data, with triplicates, RUNX2 variants indicated on the image and above the additional proteins added. Data suggest a difference in DNA binding and interaction with VDR. 1,25D3 in figure indicates 1,25(OH) 2 D 3 . doi:10.1371/journal.pone.0042617.g003 amount of GST-pull down signal in the presence or absence of 1,25(OH) 2 D 3 ligand.
Quantitative gel-shift assays were done using in vitro translated RUNX2 variants and a synthetic RUNX2 element from the mouse osteocalcin gene. These data showed a significant difference in the binding capacity of the Q-variants to the response element, suggesting that variation in the Q repeat can influence DNA binding, at least on some elements. In vitro conditions for ligand (1,25(OH) 2 D 3 ) dependent VDR binding to DNA were previously described [36]. Under these conditions, the RUNX2 response element binding activity was increased by addition of in vitro translated VDR protein and further enhanced by addition of 1,25(OH) 2 D 3 . Despite increased DNA affinity in the presence of VDR and 1,25(OH) 2 D 3 a significant difference in the in vitro DNA binding activity of Q-variants persisted such that 16Q and 30Q variants were less active than 23Q wild type (Fig. 3 C, D).

Discussion
Various lines of evidence suggest that the RUNX2 Q/A repeat region has functional polymorphism. We screened samples from various studies for Q-repeat mutations in RUNX2. A total of 26 Q-variant carriers were identified, involving 24 genetically unique individuals. The Q-variants did not present with cleidocranial dysplasia (CCD) and thus we hypothesized that the Q-variants may alter another skeletal phenotype such as BMD. Overall, Qvariants were significantly lower for BMD measures including bone ultrasound attenuation of the calcaneus, and DEXA measures of BMD at several sites. After pooling BMD measures from all studies, using Z scores, an overall estimate of the effect of Q-variant status was obtained. The genetic effect of Q-variants was around 0.7 Z-score (SD) reduction in bone density whether measured by DEXA or ultrasound. These data confirm that the Q-tract of RUNX2 has a reasonably high rate of population variation and furthermore, suggest that this variation manifests in the phenotype of lower BMD. Q-variant heterozygotes were detected at around 0.4%, or around 40,000 per ten million inhabitants. Q-variant homozygotes are expected to be infrequent; but such homozygotes may have a more severe phenotype with a mean effect of 21.4 Z-scores, if the genetic effect is additive. Whether Q-insertions have a more extreme phenotype than deletions, as suggested by the BMD data, remains to be determined in larger populations.
New repeat variants can occur through strand slippage during replication of repeat sequences such as that encoding the RUNX2 Q-repeat. Expansion of Q-repeats is known in other diseases such as Kennedy's disease where expanded androgen receptor repeats result in motor neuron disease [40] and, incidentally, the length of the repeat has been related to transcriptional function [41]. Long Q-repeats are thought to adopt b-sheet structures [42] that coil into a nanotube by forming a super helix with a critical point of increased stability at a length of 30Q [43]. This might suggest that the super structure of the Q-repeat, rather than the length, is important for altered function. Although the Q-variants observed in this study were either 30Q or less, it seems possible that longer repeat expansion of RUNX2 may be observed and may be associated with some other pathology. Interestingly, the most frequent variants observed are 16Q and 30Q, each being 7 amino acids difference from the standard wild type 23Q allele. This represents a deletion or expansion of 21 base pairs in the DNA helix, being exactly two turns of the helix increase or decrease, given the accepted value of 10.5 base pairs per helical turn in DNA.
The transcriptional activities of the 16Q and 30Q mutants were analysed using three different luciferase reporter gene assays. The variants exhibited transactivation function, however at significantly reduced levels compared to wild type (23Q). The effect was of greater magnitude using a truncated mouse osteocalcin promoter, p147, where 16Q and 30Q variant activity overlapped that of some RUNX2 mutations that are associated with cleidocranial dysplasia. The reduction in transactivation function is consistent with the decrease in BMD observed. In the case of a heterozygote with a normal wild-type allele, the net effect would be expected to reflect the activity of individual alleles. The significant effect of Qvariant RUNX2 on transactivation activity in transfection assays could result from effects on DNA binding, nuclear residence and interactions with partner proteins. Nuclear residence of Q-variants was decreased significantly in transfected cells while quantitative gel shift assays suggested a difference in binding to consensus DNA elements. We cannot distinguish between the nuclear import effects and direct effects on transcriptional machinery, since we did not examine in vitro transcription. Regardless of the mechanism, it seems plausible that the quantitative reduction in the transactivation capacity of RUNX2 is associated with overall decrease in BMD seen in human carriers. Accumulated over a lifetime, this may be enough to explain the observed effects on bone density. On the other hand, only one semi-synthetic and two natural target genes were tested, so the effect of Q-variants on other genes is not known. The vitamin D receptor (VDR) is ligand activated, is known to bind RUNX2 by interaction on the osteocalcin promoter [44], and has potent direct effects on bone cells. In the human, but not in mouse, VDR activates osteocalcin through a vitamin D response element. In this case, partial recovery of the transactivation potential of the RUNX2 Q-variants on the human osteocalcin promoter was observed on co-transfection and activation of VDR with the ligand 1,25(OH) 2 D 3 .
Are the Q-variants related to fracture? While it is suggestive that a significant effect on incident fracture was observed within the Western Australian study, this outcome should be viewed with caution due to the small numbers of Q variants. In addition, the apparent enrichment of Q-variants in post-menopausal clinic fracture cases (four variants in 303 cases genotyped) should also be considered with caution. Given the rare nature of these Qvariants, much larger studies are required to determine the relationship with fracture with precision. A frequently quoted proposition is that for every one standard deviation decrease in BMD, a person increases two fold in liability for fracture. If the effect on BMD of Q-variants is similar to the measured amount (0.7SD deficit), we can expect an increase in liability to fracture of slightly less than two-fold, if the effect of a gene such as RUNX2 on fracture acts solely via BMD. Based on the observed allele frequency, very large studies indeed (10,000 cases and controls) would be required for reasonable power test the hypothesis of association with fracture at 99% confidence for an odds ratio of around 2. If the hypothetical effect on fracture is greater than that suggested by BMD differences, then it may be possible to detect an effect in smaller numbers.
Severe RUNX2 mutations that cause CCD have a more profound effect in transfection assays, usually abolishing transactivation [6,7,19,20]. Milder phenotypes are associated with RUNX2 mutants that have some transactivation function. In some transactivation assays, notably using the p147 mouse osteocalcin promoter, the 30Q RUNX2 variant resulted in transactivation in a manner similar to CCD mutant RUNX2 forms (Q266X), suggesting a phenotype similar to CCD in terms of transfection. The subjects that carry Q-repeat mutations do not present with CCD, thus the Q-variant proteins transactivate target promoters at levels that are high enough to avoid the manifestations of CCD but at levels that reduce bone density within the normal range. Alternatively, some other mechanism may keep the phenotype out of the range of CCD, such as the CBFB or VDR interactions that were investigated in this study. The transfection data are consistent with this hypothesis. Our constructs were based on the RUNX2 protein that starts MRIPV, an isoform that is expressed off the second promoter of RUNX2. At least in mouse, and possibly in human, an isoform beginning MASNS is expressed off the first promoter. We have not tested the effect of Q-variants within the context of the isoform RUNX2-II, expressed off the first promoter. In addition, recent evidence suggests that Q-repeat regions contain an activity similar to a trans-activation domain [45,46]. Furthermore, RUNX2 Q/A repeat variation is related to skull shape in carnivores [25,26] suggesting that functional polymorphism is associated with this repeat: unfortunately no data on facial characteristics or skull shape was taken from the subjects presented in this study. It is also notable that the RUNX2 locus shows evidence of intense evolutionary selection in the comparison between Neanderthal and modern human DNA sequences [47]. How this relates to functional change remains to be established. Our data support a role for the Q-repeat in contributing to transactivation activity, although we have no data on a biochemical mechanism. It seems reasonable to conclude that a number of activities are altered to arrive at a change in target gene transactivation: including nuclear localization, DNA binding and possibly the activity of the Q-repeat as a transactivation function. The data suggests lower transcriptional activity of RUNX2, contributing to lower bone density and possibly flowon increased risk of fracture, through reduced transcriptional activity in osteoblasts.
There are several limitations of this study. No subjects were rejected from the collections for any reason prior to genotyping although some DNA failed to amplify for technical reasons. Recruitment bias was minimised as far as could be done in such clinical collections and there is no reason to suspect that genetic stratification might have influenced the outcome, although this was not formally tested. Similar results were observed in different locations, although all populations were essentially Caucasian. Bone phenotype data were derived from different commercial devices in each population and were compared through the use of standardized Z scores or Monte Carlo modelling; which should remove device-related differences. The study was limited to female Caucasians. The statistical outcome observed was not dependent on modelling covariates, although adjusted standardized values are presented. While statistically significant, these outcomes do not extend to the whole genome level of significance used in genome scans, although this study was not designed as a genome scan. The Q-variants observed in this study may be new mutations, although intergenerational studies are needed to verify this. If these are new mutants such a mutable locus may not be detected in a whole genome association study, as new mutations would be sporadically associated with genetic markers used in genome scans studies such as haplotype blocks. We have not formally established that these variants are new mutations, although the range of alleles (15Q, 16Q, 18Q and 30Q) suggests new mutations and these are compatible with strand slippage models of DNA replication through repeated sequences. Furthermore, the biochemical data based on transfections for Q-variants were tested only in 16Q and 30Q types and only within the context of RUNX2-I isoform. Biochemical effects of Q-variants may be different in the RUNX2-II isoform, which is more osteoblast specific [48] and has a differential interaction with CBFB [49] compared to RUNX2-I. Three different target constructs were tested, with similar results, although other target promoters may give different results. Finally, interactions of VDR and RUNX2 were tested in artificial constructs with a GST moiety, and tested on a single DNA element, so differences observed may not relate to the situation in the human and on other promoters.
In conclusion, the lines of evidence presented here suggest that variation in the glutamine repeat of RUNX2 has functional consequence. Significantly different phenotypic values for bone density parameters were observed in individuals heterozygous for the Q-variants. Significant functional differences were observed in assays of transcription factor activity on different target gene promoters. Furthermore, an alteration in both in vitro DNA binding activity and nuclear residence was observed, suggesting a number of different cumulative effects might occur with these variants. A partial correction of the defect in Q-variant transcriptional activity was observed in studies of interaction with two nuclear partners, CBFB and VDR, which are candidates for increasing residence of the RUNX2 variant on DNA. In addition, due to the large number of potential carriers, this genetic marker may be of some interest in the clinical setting. Intriguingly, transfection of the VDR, in the presence of active 1,25(OH) 2 D 3 overcame some of the in vivo defect in RUNX2 Q-variant, suggesting a possible route for clinical studies of Q-variant carriers or at least future studies of the bone phenotype of Q-variants relative to vitamin D status.