Identification of Candidate Gene Regions in the Rat by Co-Localization of QTLs for Bone Density, Size, Structure and Strength

Susceptibility to osteoporotic fracture is influenced by genetic factors that can be dissected by whole-genome linkage analysis in experimental animal crosses. The aim of this study was to characterize quantitative trait loci (QTLs) for biomechanical and two-dimensional dual-energy X-ray absorptiometry (DXA) phenotypes in reciprocal F2 crosses between diabetic GK and normo-glycemic F344 rat strains and to identify possible co-localization with previously reported QTLs for bone size and structure. The biomechanical measurements of rat tibia included ultimate force, stiffness and work to failure while DXA was used to characterize tibial area, bone mineral content (BMC) and areal bone mineral density (aBMD). F2 progeny (108 males, 98 females) were genotyped with 192 genome-wide markers followed by sex- and reciprocal cross-separated whole-genome QTL analyses. Significant QTLs were identified on chromosome 8 (tibial area; logarithm of odds (LOD) = 4.7 and BMC; LOD = 4.1) in males and on chromosome 1 (stiffness; LOD = 5.5) in females. No QTLs showed significant sex-specific interactions. In contrast, significant cross-specific interactions were identified on chromosome 2 (aBMD; LOD = 4.7) and chromosome 6 (BMC; LOD = 4.8) for males carrying F344mtDNA, and on chromosome 15 (ultimate force; LOD = 3.9) for males carrying GKmtDNA, confirming the effect of reciprocal cross on osteoporosis-related phenotypes. By combining identified QTLs for biomechanical-, size- and qualitative phenotypes (pQCT and 3D CT) from the same population, overlapping regions were detected on chromosomes 1, 3, 4, 6, 8 and 10. These are strong candidate regions in the search for genetic risk factors for osteoporosis.


Introduction
Osteoporosis is a common multifactorial disorder characterized by reduced bone mineral density (BMD) and compromised bone quality, with micro-architectural deterioration leading to an increased susceptibility to fracture [1]. The ability of bone to resist fracture is determined by several highly heritable phenotypes including BMD, bone size, bone structure and strength [2,3,4]. Identification of genes and pathways regulating these phenotypes and thereby underlying bone strength could provide valuable insight regarding susceptibility to osteoporosis and fracture risk.
Translational genetics starting with experimental models is a valuable tool to dissect complex genetic traits. By using inbred strains and crosses between strains, genetic heterogeneity is dramatically reduced and experimental tests, such as destructive testing of bone, can be performed. One example of successful translational genetics is the identification of the arachidonate lipoxygenase 15 (Alox15) gene encoding an enzyme that modifies polyunsaturated fatty acids, as a candidate for osteoporosis. Mice deficient in Alox15 were found to have increased bone mass [5] and subsequent association analyses in humans have shown association with BMD for SNPs in the human orthologues ALOX12 and ALOX15 [6,7].
Most experimental genetic studies on bone have been performed in the mouse [8,9,10,11,12], but observed variations in BMD, bone structure, and fragility between different inbred rat strains [13], together with progress in mapping the rat genome has established the rat as a useful genetic model for skeletal fragility [14,15,16]. Rat models offer several advantages in studying biomechanical phenotypes reflecting the ability of bone to resist fracture, since the larger bone size of rats compared to mice allows greater precision of both structural and biomechanical measurements. The biomechanical three-point bending test (3PB) used in this study provides clinically relevant phenotypes comparable to fracture in humans [17]. Fracture susceptibility is affected by a number of conditions including diabetes, where hyperinsulinemia has anabolic effects and hypoinsulinemia is associated with bone loss and increased risk of fracture [18,19,20,21,22]. Diabetes is, however, less extensively studied with regard to fracture compared to many other conditions. The inclusion of a rat strain with diabetes-associated phenotypes in linkage analysis of bone phenotypes thus offers the possibility to identify shared mechanisms e.g. by identifying overlapping quantitative trait loci (QTLs) for these traits.
In our previous studies using reciprocal F2 crosses from inbred diabetic GK and non-diabetic F344 rats, we reported results from a genome-wide screen for size, trabecular and cortical bone structure of tibia [23,24]. We identified several significant QTLs for both bone size and structure and additionally found evidence of sex-and reciprocal cross-specific interactions with bone traits. The reciprocal crosses differ with regard to grand-maternal origin, either GK or F344. This has epigenetic effects on the F2 offspring, but also results in the inheritance of different mitochondria. Mitochondria have small, circular genomes encoding 13 proteins taking part in the electron transport chain. The mitochondrial encoded proteins need to interact with proteins encoded by genes in the cell nucleus to function correctly. Polymorphisms in nuclear or mitochondrial DNA affecting such interactions could give rise to an observed reciprocal effect i.e. differing phenotype, depending on the combination of nuclear and mitochondrial alleles in the F2 population. Of note, the GK and F344 rat strains used differ in more than 100 positions in their mtDNA [25,26]. The aims of this study were to identify QTLs linked to two-dimensional dualenergy X-ray absorptiometry (DXA) and bone biomechanical phenotypes in reciprocal rat F2 crosses and to combine these with previously reported QTLs for bone size and structure. This study thus addresses all the primary skeletal determinants of fracture risk including areal BMD (aBMD), bone size, structure and strength and further delineates the involvement of sex-and reciprocal cross on the genetic regulation of bone strength related phenotypes.

Animals
Rat strains GK/KyoSwe (GK) and F344/DuCrlSwe (F344) were maintained through brother-sister breeding. Briefly, two separate F2 intercrosses were generated: one originating from grand-maternal GK and grand-paternal F344 (cross 1, (GK female6F344 male) F1), and the other from grand-maternal F344 and grand-paternal GK (cross 2, (F344 female6GK male) F1). The two groups of reciprocal F1 progeny were mated separately to yield two reciprocal F2 populations. The rats were maintained under controlled conditions as reported previously [23]. F2 progeny were sacrificed at a mean age of 215 days and body weight and length from nose tip to tail-base were recorded. The left tibia was dissected free of fat and muscle and stored in 70% ethanol at 220uC prior to biomechanical testing and analysis of skeletal phenotypes. Genomic DNA was purified from rat liver as previously described [23]. All experiments were performed with approval from the local Animal Ethics Committee; Stockholm's norra djurförsöksetiska nä mnd (N189/97), Malmö/Lunds djurförsöksetiska nä mnd (M215-01 and M143-04).

Bone strength-related phenotypes
Biomechanical testing (3PB test): Prior to testing, tibia were thawed and equilibrated at room temperature (,2 hours). Tibia were positioned on the lower supports of a three-point bending fixture (a span of ,16 mm) and held in a stable position by a 2N preload. Using an Instron 4465 materials testing machine (Instron Inc., Canton, MA, USA) with a 1KN load cell, the bones were loaded at their midpoint at a deformation rate of 1 mm/min until fracture. Load-displacement data representing structural or extrinsic properties of the bone were calculated from loaddisplacement curves and collected using LABView data acquisition software [27]. Parameters included: ultimate force (N; height of curve) reflecting the maximum load the bone can absorb before failing (i.e. bone strength), stiffness (N/mm; initial slope of the load displacement curve), and work to failure (mJ; area under curve) reflecting the total energy the bone can absorb before fracture.
Prior to mechanical testing, whole tibiae were scanned by DXA using the PIXImus TM densitometer (GE Lunar, Madison, WI). Bones were precisely positioned at the center of the X-ray cone beam and scanned in air on the Plexiglas platform provided. A grid on the Plexiglas platform allowed the position of each tibia to remain consistent between samples. aBMD was automatically calculated from the bone mineral content (BMC) and the userdefined measured region of bone area [28]. The instrument was calibrated daily using the manufacturer's phantom. Additionally, we refer to previously reported phenotypes obtained from threedimensional computed tomography (3D CT) and peripheral quantitative CT (pQCT) [23,24].

Statistical analysis
All phenotypes were normally distributed or log-transformed to obtain a normal distribution. To compare the bone phenotypes between males and females and between the reciprocal crosses, one-way ANOVA was used. The level of significance was set at p,0.05. Unless stated, p-values are nominal. Phenotypes were adjusted for reciprocal cross, age, litter size, and body-weight using regression analyses. Residuals were checked for normality and used in the QTL analysis. To account for gender attributed bone quality differences, the residuals were computed separately for each sex. Correlations between all measured bone phenotypes in the F2 progeny were evaluated using Pearson's correlation coefficients.

Genotyping and linkage analysis of quantitative traits
A total of 192 genome-wide microsatellite markers at a spacing of 10 cM, were genotyped in the F2 progeny (108 males, 98 females) and a genetic map was generated as described previously [23]. Linkage analysis was performed for each sex separately. In order to identify possible interaction differences between loci in the nuclear genome and mitochondrial DNA, the sex separated F2 progeny were also separated on the basis of reciprocal cross. QTLs on autosomes were identified employing MAP MANAGER/QTX v. b20 [29]. R/qtl was used to include the X chromosome in the linkage analysis and to confirm all reported autosomal QTLs [30]. Mapping QTLs on the X chromosome was conducted for each sex separately without further separation by cross. Permutation tests were performed to establish genome-wide significance levels by randomization of the phenotypic data in relation to genotypic data [31]. Significant (i.e. genome-wide false-positive rate of ,5%) and suggestive (i.e. genome-wide false-positive rate of ,63%) linkage was employed to establish genome-wide thresholds [32,33]. The likelihood ratio (LR) for suggestive linkage was 10.7 (logarithm of odds (LOD) = 2.3), and for significant linkage the LR range was 17.4-17.7 (LOD = 3.8). The approximate size of the QTL was defined as the region covered by a 1-LOD reduction for any of the bone traits.

Evaluation of sex-and reciprocal cross specific QTLs
Evaluation of sex-and reciprocal cross specific QTLs were performed for QTLs that reached significant linkage (LOD$3.8), following the method described previously [23]. To identify sex-specific QTLs, the LOD score differences between males and females across the genome were assessed (DLOD sex score). We applied a permutation method to evaluate sex-specific QTLs, where thresholds were established using two randomly selected equal sized subsets of males and females [34]. The randomization was conducted within each cross. Subsequently, the bone phenotypes in the two subsets were permutated to calculate DLOD sex scores across the genome. Genetic markers on the X chromosome were not included in the permutation tests. The average DLOD sex score for genome-wide significant sex-specificity at suggestive (a = 0.63) and significant (a = 0.05) level were 2.3 and 3.7, respectively.
Within each sex, subsequent reciprocal cross-separated linkage analyses were conducted to identify cross-specific QTLs. The LOD score differences between cross 1 and cross 2 (DLOD cross score) across the genome were evaluated. Thresholds of the cross specific QTLs were computed by permutation using two randomly selected equal sized cross 1 and cross 2 subsets. The average DLOD cross score for genome-wide significant cross-specificity at suggestive (a = 0.63) and significant (a = 0.05) level were 2.4 and 3.9, respectively.
To confirm the sex-and cross-specific QTLs identified with the DLOD method, likelihood ratio tests were performed comparing a full model containing a QTL6sex interaction term/cross interaction term and a reduced model without the interaction term. Both models used male and female data for sex interaction and data from the two crosses in each sex for cross interaction.
Residuals of each phenotype were examined for normality using normal probability plots. The level of significance for a specific QTL interaction with sex or cross was set at p,0.05.
A statistical power calculation for sample size, using the method of Lynch and Walsh [35] was performed as described previously [23]. Using a LOD score of 2.4 to control false positive detection of linkage, a sample size of 52 is necessary to achieve 80% statistical power for detecting a QTL with R2 value (i.e. fraction of phenotypic variance explained) of 0.25.

Influence of sex and reciprocal cross on biomechanics and DXA
Strong sexual dimorphism was observed for all biomechanical and DXA traits in the F2 progeny with significantly higher mean values in males. In contrast, no differences were observed between phenotypic mean values in the two reciprocal F2 crosses (Table 1).

QTLs linked to biomechanics
Results from the sex-and reciprocal cross separated QTL analyses for biomechanical phenotypes of tibia are summarized in Table 2. No significant linkage to biomechanical properties was identified in males when both reciprocal crosses were combined. When separating males by reciprocal cross, a QTL on chromosome 15 linked to ultimate force reached genome-wide significance in cross 1 with GK grand-maternal origin. This QTL met the criteria for genome-wide cross-specificity at the suggestive level (DLOD cross $2.4), and the likelihood ratio test confirmed a cross-specific interaction for this locus (LR = 10.4, p = 0.005) (Fig. 1B). An overlapping QTL in females from cross 2 was suggestively linked to work to failure ( Table 2).
The maximum LOD-score for biomechanical phenotypes in females was 5.5 on chromosome 1 (18-43 cM) with significant linkage to stiffness. This locus also showed linkage to ultimate force at the suggestive level. The chromosome 1 QTL linked to stiffness reached genome-wide significance for female-specific interaction with DLOD sex evaluation (DLOD sex $3.7, results not shown) but did not reach significance in the likelihood ratio test. Suggestive QTLs for biomechanical phenotypes are reported in table S1.

QTLs linked to DXA phenotypes
Results from the sex-and reciprocal cross separated QTL analysis for DXA traits are summarized in Table 2. In all males, significant linkage to BMC and to area was identified on chromosome 8 (35-59 cM). The chromosome 8 BMC QTL indicated male-specific   interaction (DLOD sex .2.3, results not shown) but was not confirmed in the likelihood ratio test. In males with F344 grand maternal origin (cross 2), significant linkage was detected on chromosome 2 for aBMD, and on chromosome 6 for BMC. The QTL on chromosome 2 also overlapped with a suggestive QTL linked to BMC in cross 2 males. The significant QTLs on chromosome 2 (aBMD) (Fig. 1A) and 6 (BMC) reached genome-wide significance in the DLOD cross evaluation (DLOD cross $3.9) and cross-specific interactions were confirmed for both by likelihood ratio tests (LR = 10.4, p = 0.005 and LR = 10.9, p = 0.004) respectively. In females, QTL analysis on the X chromosome without separation by cross identified a significant QTL linked to tibial area. The cross-separated QTL analysis revealed significant linkage to BMC on chromosome 4 (36-48 cM) in females with GK grand maternal origin (cross 1). An overlapping region also showed suggestive linkage to aBMD in cross 1. The QTL identified on chromosome 4 significantly linked to BMC reached genome-wide significance in the DLOD cross evaluation (DLODcross $3.9), but did not reach significance for cross-specific interaction in the likelihood ratio test. Suggestive QTLs for DXA phenotypes are reported in table S1.

Co-allelic effects
The GK allele had an increasing effect on genotypic mean values for the majority of the QTLs significantly linked to both biomechanical and DXA traits in males and females (Table 3). There were however also heterosis-like effects, where the heterozygous group carrying one GK-and one F344 allele differed significantly from both homozygous groups and had higher mean aBMD (marker D2Mgh5) and ultimate force (marker D15Mit2).

Co-localization of QTLs for bone strength-related phenotypes
By combining the QTLs identified in this study with QTLs linked to other bone-related phenotypes, overlapping or colocalized QTLs for biomechanical properties (3PB), bone density (DXA) and/or bone quality (pQCT and 3D CT) phenotypes were detected on chromosomes 1, 3, 4, 6, 8 and 10 (Table 4). In order to increase understanding of the observed genetic co-localizations, correlation coefficients between phenotypes mapping to the same region were evaluated. A large region on chromosome 1 (17-79 cM) displayed linkage to multiple phenotypes obtained from all four methods (3PB, DXA, pQCT and 3D CT) in both males and females. The biomechanical phenotypes ultimate force and stiffness were strongly correlated (0.64-0.88) with the other phenotypes linked to this region except for cortical volumetric BMD (cortvBMD) and metaphyseal volumetric BMD (me-tavBMD;20.02-0.24). Similar patterns were seen at other loci, where most co-mapped phenotypes were correlated with each other, with exceptions mainly attributed to biomechanical phenotypes. For phenotypes linked to chromosome 3 (3-24 cM), ultimate force strongly correlated with 3D CT and pQCT phenotypes (0.71-0.88), while correlations with total bone mineral content (totBMC) and metaphyseal volumetric bone mineral density (metavBMD) were low (0.04-0.14). The QTLs that colocalized to chromosome 4 (36-74 cM) were measured by 3PB, 2D DXA and 3D CT and these phenotypes were strongly correlated (0.63-0.96) with the exception of work to failure which demonstrated low correlation with the other phenotypes (0.07-0.26). The same pattern was observed for chromosome 8 (35-63 cM) with overlapping QTLs from all four methods, but moderate correlations between work to failure and the other comapped phenotypes. In contrast, strong correlations (0.79-0.91) were observed for all phenotypes that co-localized to chromosome 6 (25-60 cM) and to chromosome 10 (73-91 cM) (0.45-0.93), regions where biomechanical phenotypes did not co-map.

Discussion
In this study, we identified multiple QTLs for biomechanical and two-dimensional DXA phenotypes in F2 progeny of GK and F344 rats and confirmed bone QTL interactions with reciprocal cross. Previously, we reported results from a genome-wide screen for bone size and trabecular and cortical bone structure of tibia in the same F2 progeny [23,24]. The biomechanical testing reported here provides bone phenotypes directly reflecting the ability of bone to resist fracture. We have thus addressed the majority of clinically important bone properties related to bone strength and fracture susceptibility.
By combining QTLs linked to distinct bone phenotypes measured by separate methods, overlapping chromosomal regions linked to bone size, structure and strength were identified. Two large regions on chromosomes 1 (17-79 cM) and 8 (35-63 cM) displayed linkage to multiple phenotypes obtained from four different methods (3PB, DXA, pQCT and 3D CT) in both males and females. The regions with overlapping QTLs influenced different phenotypes, the majority of which were correlated at a   population level. Exceptions mostly represented biomechanical phenotypes. As an example, overlapping QTLs on chromosome 3 included the biomechanical phenotype ultimate force which was strongly correlated with the dimensional phenotypes 3D CT-and pQCT but only weakly correlated with totBMC and metavBMD. The co-mapping of these phenotypes was thus not due to phenotype correlation, but represents different aspects of the QTLs. Co-localization of QTLs for uncorrelated phenotypes could reflect the existence of distinct sub-loci. Conversely, highly correlated co-mapped bone phenotypes could reflect pleiotropic gene effects with dependence between phenotypes. Thus, the observed co-localizations of QTLs are unlikely to be explained by phenotypic correlation alone and the regions could each contain either a single QTL with pleiotropic effects, or several QTLs that influence the correlated phenotypes independently. Combining QTLs from different methods allows greater characterization of the gene regions, indicating whether they are likely to contain one or several candidate genes and if these are likely to participate in the same biological processes. Since the biomechanical properties of bone are the sum of many different characteristics, the lower correlation of these traits to specific quantitative and qualitative traits was expected. The notion that QTLs for bone phenotypes interact specifically with sex have been reported previously [8] but the observed interactions between nuclear QTLs for several bone phenotypes and reciprocal crosses, that differ with regard to mitochondrial DNA sequence, demonstrates a new and important aspect to be considered when interpreting the genetics of phenotypes related to bone strength.
In this study, the reciprocal cross-separated QTL analysis in males allowed identification of cross-specific interactions for three QTLs on chromosomes 2 (aBMD), 6 (BMC) and 15 (ultimate force). The QTLs on chromosomes 2 and 6 were identified in all males at a suggestive level but showed genome-wide significant linkage only in the cross with F344 grand maternal origin, while the QTL on chromosome 15 was only detected in the cross with GK grand maternal origin and would not have been detected in the combined sample including both reciprocal crosses. This illustrates the importance of sub-group analysis and of considering reciprocal cross effects in order to improve the detection of candidate genes for fracture susceptibility.
The respective reciprocal cross did not affect the phenotypes at a population level, as no significant differences were seen between the phenotypic mean values between cross 1 and cross 2 in either males or females (Table 1). Instead, the reciprocal effects were seen as interactions with specific QTLs and thus provide support for nuclear-mitochondrial interactions to be involved in the genetic regulation of bone phenotypes. Notwithstanding the mounting evidence for mitochondrial interactions involved in complex diseases such as osteoporosis and type-2 diabetes, other factors such as genomic imprinting, the presence of QTLs on sex chromosomes and maternal environment, are theoretical explanations for reciprocal cross specific inheritance and may contribute to the observed effects as previously discussed [23]. Additionally, since all cross-specific QTLs in this study were identified in males, we cannot exclude the possibility that interactions with a Y chromosome linked locus could explain the observed reciprocal cross effect. However, the more than 100 variant positions in both coding and non-coding regions that have been identified in GK mitochondrial DNA compared to F344 [25,26], support mitochondrial genotype as a possible factor behind the observed reciprocal cross effect in this study. The reciprocal cross interactions could thus reflect functional interaction between nuclear-and mitochondrial encoded proteins.
Although the candidate regions are large and contain many genes, it is interesting to note that the QTL regions on chromosome 4 and 6 that demonstrated significant reciprocal cross interaction contain genes encoding mitochondrial proteins. These include mitochondrial ribosomal protein S33 [Mrps33], the NADH dehydrogenase (ubiquinone) 1 beta subcomplex 2 [Ndufb2] and glutathione S-transferase kappa 1 [Gstk1] on chromosome 4 and ATP synthase subunit s [Atp5s] and Acylcoenzyme A thiosterase 2 [Acot2] on chromosome 6. In addition, the region on chromosome 4 overlaps a previously identified QTL strongly linked to femoral neck structure phenotypes in female (F3446LEW) F2 rats [14].
Several QTLs identified in this study co-localized with our previously reported QTLs linked to tibial bone phenotypes obtained from pQCT and 3D CT [23,24] and were in many cases influenced by sex-and reciprocal cross ( Table 4). The region on chromosome 1 (17-79 cM) linked to multiple phenotypes obtained from all four methods (pQCT, 3D CT, 3PB and DXA) in both males and females showed predominantly higher significance in the reciprocal cross 1 carrying GK mtDNA. Several QTLs linked to cortical bone traits in other combinations of rat strains have been mapped to this region [14,16], emphasizing the importance of this large chromosome region for variation in bone traits between individuals. The osteoporosis candidate genes transforming growth beta 1 (TGFB1) [36] and estrogen receptor alpha (ESR1) [37] are both localized within this locus. Interestingly, QTLs for fasting glucose also map to this region [38] making it a strong candidate for focused investigation of possible shared or linked polymorphisms regulating diabetes-and osteoporosisrelated phenotypes. To clarify shared mechanisms between diabetes and osteoporosis, such analyses could include congenic strains and/or advanced intercross lines. It is recognized, as a possible limitations, that this chromosomal region is large and likely to harbor several sub-loci.
Several established osteoporosis candidate genes are also found within the human regions syntenic to the chromosome 4 QTL linked to multiple rat bone strength phenotypes. These include the calcitonin receptor (7q21.3), collagen 1 alpha 2 (7q21.3) and WNT 2 and WNT 12 from the WNT signalling pathway. The cytokine macrophage migration inhibitory factor (MIF) gene, recently shown to be associated with bone loss in elderly women [39] is also located within this region. The human regions syntenic to the rat chromosome 6 QTL contains around 140 genes. This region is too large to pinpoint candidate genes, but interestingly contains genes which contribute to skeletal development e.g. the transcriptional regulator TWIST1 (7p21) and GDF7 (2p24), a member of the TGFB superfamily. Further delineation of the complex genetic architecture of bone phenotypes for fracture susceptibility would ideally involve a genome screen for epistatic interactions between the identified loci. However, such analysis requires a larger sample size than was available in this study. Congenic strains or advanced intercross lines are possible tools for future fine mapping of the identified QTLs, and would allow for positional cloning of candidate genes and epistatic interaction analyses, respectively.
In summary, we have identified QTLs for biomechanical and twodimensional DXA phenotypes influencing bone strength with interaction from reciprocal cross in an F2 intercross between GK and F344 rats. The observed interaction between nuclear QTLs and reciprocal cross for numerous bone associated phenotypes supports the potential for mitochondrial effects on bone. By combining data from this and our previous studies, we were able to identify specific but also overlapping chromosomal regions for bone size, structure and strength. These findings illustrate the importance of analysing different determinants contributing to bone strength in order to detect candidate genes or pathways for bone regulation, from the macro-to the micro-structural level. Of particular interest is the identified QTL region on chromosome 1 linked to many bone phenotypes and also reported to affect fasting glucose. This region is thus a strong candidate for the identification of genes contributing to bone regulation and potentially type-2 diabetes. Table S1. QTL size defined as the region covered by a 1-LOD reduction for any of the bone traits. a Genome-wide suggestive QTLs (LOD$2.3) are marked in bold. (DOC)