Genetics of evolved load resistance in the skeletons of unusually large mice from Gough Island

Abstract A primary function of the skeleton is to resist the loads imparted by body weight. Genetic analyses have identified genomic regions that contribute to differences in skeletal load resistance between laboratory strains of mice, but these studies are usually restricted to 1 or 2 bones and leave open the question of how load resistance evolves in natural populations. To address these challenges, we examined the genetics of bone structure using the largest wild house mice on record, which live on Gough Island (GI). We measured structural traits connected to load resistance in the femur, tibia, scapula, humerus, radius, ulna, and mandible of GI mice, a smaller-bodied reference strain from the mainland, and 760 of their F2s. GI mice have bone geometries indicative of greater load resistance abilities but show no increase in bone mineral density compared to the mainland strain. Across traits and bones, we identified a total of 153 quantitative trait loci (QTL) that span all but one of the autosomes. The breadth of QTL detection ranges from a single bone to all 7 bones. Additive effects of QTL are modest. QTL for bone structure show limited overlap with QTL for bone length and width and QTL for body weight mapped in the same cross, suggesting a distinct genetic architecture for load resistance. Our findings provide a rare genetic portrait of the evolution of load resistance in a natural population with extreme body size.


Introduction
The skeleton supports the body, protects vital organs, and facilitates locomotor movement.Bone stiffness and strength are crucial to these functions.Bones must resist deformation (stiffness) and failure (strength) to handle the loads caused by bearing weight and performing movement (Jepsen 2009).
The ability of bone to withstand mechanical loads is determined by its material and structural properties (Hart et al. 2017).On the material side, the quality of bone is primarily influenced by its mineral density, which is the product of mineralization and porosity (Hart et al. 2017).From the structural perspective, cross-sectional geometry is a key determinant of a bone's ability to resist loads (Adams and Ackert-Bicknell 2015).The most relevant morphological parameter for loading in compression (pushing) or tension (pulling) is a cross-sectional cortical bone area, whereas the cross-sectional moment of inertia describes the potential for resistance to bending and torsion (Jepsen et al. 2015).
Inbred strains of mice raised in the same laboratory environment differ in structural and material properties of long bones (Beamer et al. 1996;Jepsen et al. 2003;Wergedal et al. 2005;Reich et al. 2008), demonstrating the contribution of genetic factors and positioning mice as a powerful experimental system for characterizing them.Suites of quantitative trait loci (QTL) shape variation among mouse strains in bone strength, cross-sectional geometry, and mineral density (Beamer et al. 2001(Beamer et al. , 2012;;Klein et al. 2002;Adams and Ackert-Bicknell 2015;Lu et al. 2019).QTL for bone geometry and material composition sometimes colocalize with QTL for whole-bone mechanical strength evaluated in functional tests (Jepsen 2009;Adams and Ackert-Bicknell 2015) and with QTL for obesity (Reich et al. 2008), providing genetic confirmation of the relevance of bone structure for load resistance.Fine-mapping and systems genetic analyses have identified candidate genes for traits involved in load resistance, including Qsox1, which controls cortical bone mass and strength in diversity-outbred mice (Al-Barghouthi et al. 2021).
Despite this progress, 2 important questions remain largely uninvestigated in the genetic characterization of load resistance abilities.First, are the structures that determine load resistance abilities in different bones generated by the same set of loci?Although mechanical and developmental considerations suggest potential for distinct genetic architectures across the skeleton (Hall 2005), most genetic mapping studies of bone strength focus on the femur (or the femur and the tibia), leaving this question unaddressed.Second, what is the genetic basis of evolved differences between natural populations in bone structures that resist load?Genetic mapping studies to date are restricted to classical inbred mouse strains and their derivatives, leaving open the question of how bone structure adapts to increases in load in response to novel environments in the wild.
Natural populations that have recently evolved substantial increases in body size are promising targets for the genetic investigation of skeletal load resistance abilities.Mice on Gough Island (hereafter "GI"), a remote island located near the middle of the South Atlantic Ocean, are the largest wild house mice on record (Rowe-Rowe and Crafford 1992; Gray et al. 2015), weighing close to twice as much as their mainland relatives (Jones et al. 2003;Gray et al. 2015).This remarkable expansion in size likely happened during the few hundred generations following the colonization of the island (Rowe-Rowe and Crafford 1992).GI mice stand as a prominent example of the island rule, a taxonomically broad pattern among vertebrates in which island populations evolve extreme body sizes compared to their mainland counterparts (Foster 1964;Van 1973;Lomolino 1985;Benítez-López et al. 2021).The bones of GI mice must resist the additional load imposed by their extra body weight, which could have led to the evolution of bone structural and/or material properties in these mice.GI mice belong to the same subspecies as classical inbred strains (the Western European house mouse, Mus musculus domesticus) (Gray et al. 2014), facilitating the genetic examination of skeletal traits.In previous work, we uncovered the genetic basis of evolution in body size and skeletal shape in GI mice (Gray et al. 2015;Parmenter et al. 2016Parmenter et al. , 2022)).
In this study, we use GI mice as a model system to understand how morphological determinants of skeletal load resistance abilities evolve in nature.We identify QTL that contribute to differences between GI mice and a mainland reference strain in the structures of 7 bones.Our findings extend genetic understanding of load resistance across the skeleton and provide clues into the changes that accompany the evolution of extreme body size.

Mice
GI is part of the United Kingdom Overseas Territory of Tristan da Cunha located in the South Atlantic Ocean, approximately halfway between South America and South Africa (40° 19′S and 9° 55′W).Fifty mice live-trapped on GI in September 2009 were transferred to the Charmany Instructional Facility in the School of Veterinary Medicine at the University of Wisconsin-Madison.Upon their arrival, 46 mice (25 female and 21 male) were used to establish a breeding colony.
All mice in this study were housed at the University of Wisconsin-Madison Charmany Instructional Facility (Madison, WI).Female mice and male mice were housed separately in microisolator cages with a maximum of 4 mice per cage.Ground corn cobs (1/8th inch; Waldschmidt and Sons, Madison, WI) were used as bedding; nesting material and irradiated sunflower seeds (Harlan Laboratories, Madison, WI) were provided for enrichment.The room was temperature controlled (68-72°F) and set on a 12-hour light/dark cycle.Water and rodent chow (Teklad Global 6% fat mouse/rat diet; Harlan Laboratories, Madison, WI) were provided ad libitum.Breeding individuals were fed breeder chow (Teklad Global 19% protein/9% fat; Harlan Laboratories, Madison, WI) ad libitum.All mice were weaned between 3 and 4 weeks of age.
Two partially inbred lines of GI mice were generated through full-sib mating for 4 filial generations (Gray et al. 2015;Parmenter et al. 2016).Mice from the fully inbred, wild-derived strain WSB/EiJ (subsequently "WSB") were ordered from Jackson Laboratories (Bar Harbor, ME), mated to generate mice born in our facility, and treated as a mainland reference for comparisons with GI mice.
To characterize phenotypic differences between GI mice and WSB mice, we separately constructed a linear model for each trait, with strain and sex treated as fixed effects.Sample sizes for these comparisons were 30 GI mice (17 females; 13 males) and 24 WSB mice (13 females; 11 males).
One pair of female and male siblings from each of 2 partially inbred lines of GI mice was crossed with WSB.F1s were intercrossed to generate F2s.Altogether, 4 independent F2 intercrosses were performed (Gray et al. 2015;Parmenter et al. 2016).Genetic analyses of bone traits focused on 760 F2 mice: 405 from Cross A (WSB × GI = 243; GI × WSB = 162) and 355 from Cross B (WSB × GI = 203 and GI × WSB = 152).All F2 mice were weighed at 16 weeks of age and euthanized by CO 2 asphyxiation or by decapitation.This study was approved by the Institutional Animal Care and Use Committee at the University of Wisconsin-Madison.

Phenotyping
Our study focused on limb bones including the femur, tibia, scapula, humerus, radius, and ulna, as well as the mandible.We did not survey the pelvis because identifying a cross-sectional middle is more challenging.After dissecting away soft tissues, skeletal elements were scanned by micro-computed tomography (micro-CT) at a 20.5 μm voxel resolution in a vivaCT 70 scanner (Scanco) at 70 kVp/114 mA.We output scans as stacked 2,048 × 2,048 tiff images.These images were cropped for each skeletal element in Adobe Photoshop.We imported stacked tiff images into Avizo where digital reconstructions of skeletal elements were resliced at the diaphyseal midpoint for long bones, the scapular neck, and the mandibular M1 for cross-sectional analysis.Resliced tiff images were imported into ImageJ2 (Fiji) Lines indicate orientations for calculating maximum area moment of inertia (Imax) and minimum area moment of inertia (Imin).Polar moment (J) is computed as Imax + Imin.(Rueden et al. 2017) for cross-sectional analysis using BoneJ (Doube et al. 2010).
We characterized cross-sectional properties of bones related to load resistance abilities using 4 standard measurements (Fig. 1).Cortical area (CA), an indicator of the ability to resist shear, was measured as the total area of bone in a cross-section.Imax, the maximum area moment of inertia and an indicator of maximum resistance to bending, was measured as the sum of bone area times its squared distance from the bone centroid.Imin, the minimum area moment of inertia and an indicator of minimum resistance to bending, was measured similarly and is perpendicular to Imax.J, the polar moment of inertia and an indicator of resistance to torsion, was measured as the sum of Imax and Imin.Our genetic analyses considered CA, Imax, Imin, and J collected from each of the 7 bones analyzed.
Since differences in bone quality between GI mice and WSB mice would impact the interpretation of bone cross-sectional data, we compared bone density in these 2 groups.Bone density was calculated by converting grayscale values from bone crosssections to mg/HA based on a linear regression equation created from known calibrations between grayscale values and mg/HA using manufacturer-supplied phantoms (R = 0.999): mgHA/cm 3 = −222.813+ 6.927 (grayscale value).

Genotyping
Genotypes for QTL mapping were taken from Gray et al. (2015).Briefly, livers from F2s, GI parents of the cross, and WSB parents of the cross were sent to Geneseek (NeoGene Corporation) for DNA extraction and genotyping.Mice were genotyped using the Mega Mouse Universal Genotyping Array (MegaMUGA; Geneseek, Lincoln, NE), an Illumina array platform containing 77,800 single nucleotide polymorphisms (SNPs), evenly spaced throughout the genome (Threadgill and Churchill 2012).After data cleaning, analyses focused on 11,833 fully informative SNPs that were fixed in the 4 GI mouse parents and therefore segregated as in a standard F2 intercross between inbred lines.We estimated genetic distances between SNPs assuming a genotyping error rate of 0.2% and the Carter-Falconer mapping function (Carter and Falconer 1951).The genetic map was reconstructed using 1,212 F2s (Gray et al. 2015).

QTL mapping and characterization
To identify QTL for each trait separately, we first conducted a single-QTL scan, which could detect up to one QTL per chromosome, using Haley-Knott regression in the scanone function from the R/qtl package (Broman et al. 2003).A QTL was deemed significant if its logarithm of odds (LOD) score exceeded the genomewide 5% threshold set by 1,000 permutations (Churchill and Doerge 1994).Next, we performed a multiple-QTL scan with a stepwise search using a penalized LOD score as the criterion (Broman and Speed 2002), implemented by stepwise in R/qtl.We initialized the search with QTL detected by single-QTL mapping, used the significance threshold from single-QTL mapping as the penalty, and assumed QTL effects combine additively across loci.Based on results from single-QTL scans, we set the maximum number of possible QTL at 10 for all traits except for the humerus, for which we set the maximum number of QTL at 15. Additive effects of QTL (a) were calculated as half the difference between genotypic means of GI homozygotes and WSB homozygotes and scaled by dividing by the F2 phenotypic standard deviation for the trait of interest.Dominance effects (d ) were calculated as the difference between the genotypic mean of GI/WSB heterozygotes and the midpoint of the genotypic means of the 2 homozygotes and scaled by dividing by the additive effect (d/a).
Standard errors for d/a were computed using the delta method (Lynch and Walsh 1998).The total percent of F2 variance explained by QTL was calculated from the full model.The percent of F2 variance explained by individual QTL was computed by comparing the full model to a sub-model without that QTL.QTL effects were estimated using fitqtl in R/qtl.
To map QTL while accounting for correlations among traits, we conducted principal component analyses (PCAs) of phenotypes across F2s.For these analyses, we excluded all measurements from the mandible and Imin measurements from all bones to reduce redundancy and focus on the postcranial skeleton, leaving 18 measurements from across 6 bones.We computed principal components from the pairwise correlation matrix of these 18 traits using princomp in R. We separately used scores for principal components 1 and 2 as traits in single-QTL and multiple-QTL scans.
All QTL analyses included sex and mother as additive covariates.After preliminary single-QTL scans revealed limited evidence for the effects of the X chromosome, QTL analyses were restricted to the autosomes.We estimated 1.5 LOD drop intervals for all QTL using lod.int in R/qtl.To estimate genomic positions of QTL, we found the locations of the closest genotyped SNPs in the grcm39 version of the mouse reference genome (using updates kindly provided by Karl Broman: https://raw.githubusercontent.com/kbroman/MUGAarrays/master/UWisc/mm_uwisc_v2.csv).
For a subset of QTL with overlapping confidence intervals across many traits, we compared the fit of a null model assuming a single QTL affects the set of traits to the fit of models assuming 2 linked QTL affect the set of traits.This analysis was implemented in the qtlpvl package in R/qtl using the "testpleio.1vs2"command with the "complete" search method and 100 parametric bootstrap simulations to estimate P-values under the null model of a single QTL.
We searched for candidate genes for QTL that affect many load resistance traits.First, we re-estimated the QTL position by examining all traits associated with the QTL using the "scanone.mvn"command in the qtlpvl package.In Mouse Genome Informatics (Ringwald et al. 2022), we searched the 1.5-LOD interval for genes associated with "abnormal limb long bone morphology" (MP: 0011504) and examined the resulting references for reports of bone phenotypes.We also searched the interval for genes with nonsynonymous differences between GI mice and WSB mice used in the cross.For this purpose, we applied Liftover in the UCSC Genome Browser (http://genome.ucsc.edu) to locate the 1.5-LOD interval in the grcm38 assembly.We used custom tracks in the University of California Santa Cruz (UCSC) Genome Browser we created from genome sequences of the partially inbred GI mice and WSB mice (Nolte et al. 2020) used as parents of the cross and the Variant Annotation Integrator (UCSC) to identify all nonsynonymous differences in the QTL interval that are fixed between these GI mice and WSB mice.We used the Variant Table in Ensembl to predict the deleteriousness of these nonsynonymous changes with the appproach called Sorting Intolerant from Tolerant (SIFT; https://ensembl.org/info/genome/variation/prediction/protein_function.html).

Higher skeletal load resistance in GI mice compared to mainland reference strain
Analysis of variance (ANOVA) treating strain and sex as fixed effects factors indicates that most structural traits have significantly higher values in GI mice vs WSB mice (Table 1).Relative differences are generally higher for Imax, Imin, and J than for CA.Differences vary among bones, with traits from the scapula, humerus, and mandible tending to show higher relative values.ANOVA indicates that the sexes do not differ significantly for most structural traits in GI mice and WSB mice (Table 1).One exception is the mandible, where measurements tend to show higher values in females (Table 1).
ANOVA also demonstrates that strains do not differ significantly in the density of any of the 7 measured bones (P > 0.05 for all bones).Females and males differ significantly in the density of the mandible (P = 0.007), scapula (P = 0.004), radius (P = 0.048), and humerus (P = 0.035), with females having higher mean densities for each bone.

Phenotypic correlations among F2s
Pairwise trait correlations across F2s show several patterns (Fig. 2).Pearson's correlations are universally positive with low P-values (all P ≤ 0.004).Different traits within bones are usually more highly correlated than the same measurements across bones.Femur with tibia and femur with humerus show the strongest correlations among bones.Correlations between the scapula (Imax and J) and femur, tibia, and humerus are also high.Overall, J traits have the highest correlations, followed by Imax, CA, and Imin.The weakest correlations involve the mandible or radius Imin.
To gain preliminary insights into multivariate patterns among the F2s, we conducted PCA on the pairwise correlation matrix formed from 18 traits (see Materials and methods).Principal component 1 (PC1) explains 72% of the phenotypic variance among F2s.The 18 traits included in the PCA make similar contributions to PC1, with loadings ranging from 0.21 to 0.25 across traits.PC1 score is highly correlated with 16-week body weight (Fig. 3; Pearson's r = 0.78; P < 2.2e−16).Pairwise comparisons with weights taken at other ages (Gray et al. 2015) show that the high correlation between PC1 and weight is visible earlier in ontogeny (e.g.r = 0.72; P < 2.2e−16; 16-week PC1 vs 6-week body weight).
PC2 explains 9% of the phenotypic variance among F2s.Loadings along PC2 vary among traits in both sign and magnitude.The 3 measurements from the ulna have the strongest effects; traits from the scapula and humerus have weak effects.There is a contrast in signs among groups of bones, with increases in the femur and tibia increasing the PC2 score and increases in the ulna and radius decreasing the PC2 score.PC2 score shows a much weaker correlation with 16-week body weight than does PC1 score (Fig. 3; r = 0.17; P < 4.4e−06).

QTL for skeletal load resistance traits
Single-QTL mapping and multiple-QTL mapping identify QTL for every trait except radius Imin (Supplementary Fig. 1; Table 2; Fig. 4).QTL are detected on every autosome except chromosome 17.The number of detected QTL ranges across phenotypes from 1 (Mandible CA) to 13 (Humerus Imax and Humerus J).Although we conducted no formal tests designed to gauge the extent of shared QTL effects across phenotypes, comparison of QTL positions among traits suggests variation in the breadth of QTL detection (Fig. 4; Supplementary Fig. 1).Some QTL appear to be specific to certain bones, such as QTL located on chromosomes 5 and 19 that are only detected using data from the humerus.Some QTL are detected globally, such as a QTL located on chromosome 10 that impacts every bone surveyed.Most QTLs are detected for a subset of bones.QTL often affect multiple traits, but not always, raising the possibility that resistance measures sometimes have independent genetic bases within the same bone.
Multiple-QTL mapping identifies 6 QTL that contribute to the PC1 score (Table 2).These QTL are a subset of those found across the collection of traits (Fig. 4).Although QTL for PC1 score map to some of the  (Gray et al. 2015), the 1.5-LOD confidence intervals overlap only for the QTL on chromosome 7.This finding suggests that cross-sectional traits and body weight have distinct genetic architectures, despite their functional relationship.Multiple-QTL mapping identifies 3 QTL that contribute to the PC2 score, one of which overlaps with a QTL for the PC1 score.As a group, QTL for PC1 and QTL for PC2 have narrower 1.5-LOD intervals than most measurements mapped on the original scale.
To further examine the breadth of QTL effects, we compared the fit of a null model assuming a single QTL affects a set of traits to the fit of models assuming 2 linked QTL affect the same set of traits.We focused on 3 QTL with overlapping confidence intervals across many traits: a region on chromosome 7 linked to 18 traits, a region on chromosome 8 linked to 20 traits, and a region on chromosome 10 linked to 17 traits.For the regions on chromosomes 7 and 8, we failed to reject the null hypothesis of a single QTL (P > 0.05).For the region on chromosome 10, we rejected the null hypothesis of a single QTL (P = 0.01).Two QTL separated by 5.9 cM are detected, with the most proximal QTL (at 11.35 cM) affecting 14 postcranial traits, and the other QTL (at 17.24 cM) modulating mandible Imax, mandible Imin, and mandible J.These results support the notion that some QTL impact resistance abilities across many dimensions and bones, though the power to reject the null hypothesis in favor of multiple closely linked QTL is constrained by the sample size and the number of crossovers in an F2.Standardized additive effects of QTL are summarized in Fig. 5 and Table 2. Across QTL (excluding those mapped for principal component scores and those with absolute scaled dominance effects (|d/a|) > 1), the absolute values of additive effects range from 0.149 to 0.433 F2 phenotypic standard deviations, with a mean of 0.223.Among these QTL, 93 have positive additive effects, indicating that the GI mouse allele increases the trait value, whereas 37 have negative additive effects, indicating that the GI mouse allele decreases the trait value.Each bone features QTL with a mixture of positive and negative additive effects.ANOVA treating bone and trait (CA, Imax, Imin, and J) as fixed effect factors (again excluding QTL with |d/a| > 1) finds evidence for differences in standardized additive effect absolute values among bones (F 6,120 = 2.93; P = 0.01) and among traits (F 3,120 = 2.73; P = 0.05).A Tukey post hoc test indicates that tibia-radius (P = 0.02), ulna-tibia (P = 0.03), and Imin-Imax (P = 0.03) contrasts contribute to the significant differences in additive effects among bones and among traits.
The percent variance among F2s explained by individual QTL ranges from 1.13 to 8.94% (Table 2), with a mean of 2.40%.

Candidate genes for skeletal load resistance
We focused our search for candidate genes on the shared QTL on chromosomes 10 and 7, which affect 14 postcranial traits and 18 traits (respectively) and have 1.5-LOD intervals that span less than 10 Mb (Table 3).Ccn2 is the only gene located in the chromosome 10 shared QTL interval that is known to play a specific role in the formation of the skeleton.Phenotypes of knockout mice demonstrate that Ccn2 (also known as Ctgf) regulates extracellular matrix remodeling that affects endochondral bone ossification (Ivkovic et al. 2003) (though it should be noted that long bone in adults is mostly produced by intramembranous ossification).Ccn2 harbors a single nonsynonymous SNP between the parents of our cross.The variant is located at amino acid position 1 in one transcript and amino acid position 24 in a second transcript.In one transcript, the GI variant constitutes an alternative start codon that encodes threonine instead of methionine.This substitution is predicted to be deleterious (SIFT score = 0.01).
The shared QTL region on chromosome 7 contains several genes with established effects on skeletal morphology (Table 3).Scube2 mouse mutants show reduced osteoblast cell number (Lin et al. 2015).Scube2 contains 2 nonsynonymous SNPs between the parents of our cross.Both SNPs have SIFT scores of 1, predicting no deleterious consequences.Knockout mice for Spon1 exhibit increased bone mass in the femur and tibia (Palmer et al. 2014).The single nonsynonymous SNP between the parents of our cross found in Spon1 is predicted to be tolerated (SIFT score = 0.9).Long bones from Pth knockout mice have higher cortical thickness (Miao et al. 2002).Xylt1 regulates early chondrocyte maturation in mice (Mis et al. 2014).

Discussion
Our findings point to several conclusions about the evolution and genetics of determinants of bone strength.GI mice  among traits in F2s suggest particularly tight integration between changes to the femur and changes to the tibia/humerus, which could reflect shared developmental trajectories between these bones.Although it is difficult to distinguish between pleiotropy and linkage as explanations for the colocalization of QTL in an F2 intercross, overlapping confidence intervals and mapping of principal component scores imply that several QTL alter the structure of multiple bones.Loci on chromosomes 7, 8, and 10 affect a plethora of traits throughout the skeleton.At the same time, structural QTL vary in their breadth of detection.Some QTL appear restricted to single bones and some QTL were detected for certain traits but not others within the same bone.We caution that a full picture of the distribution of shared QTL effects across bones and traits awaits analyses that consider QTL effects in a multivariate context.Nevertheless, these findings significantly expand the view of the genetics of load resistance in the mouse skeleton beyond previous investigations, which focus primarily on the femur.Our results underscore the need for genetic studies to consider multiple bones.
We detected no QTL that alter bone structure in major ways, despite the precipitous evolution of body size in GI mice.All QTL additive effects are smaller than 0.44 F2 standard deviations and each QTL explains less than 9% of F2 variance, despite an expected upward bias in effect size estimates for detected QTL (Beavis 1994).Therefore, the path to increased load resistance abilities in GI mice followed expectations for typical quantitative traits, involving many mutations with individually modest effects.A similar genetic inference has been drawn for the evolution of other morphological traits in GI mice (Gray et al. 2015;Parmenter et al. 2016Parmenter et al. , 2022)).Parmenter et al. (2016) report QTL for bone length and width in the same cross we analyzed here.Most QTL from the 2 studies do not overlap, indicating that the evolution of cortical structure and the evolution of whole-bone measures were mostly accomplished by different genes.A subset of bone structure QTL overlap with those from Parmenter et al. (2016), including QTL on chromosome 7 for humerus CA, humerus Imax, humerus Imin, and humerus J (humerus midshaft diameter) and for tibia CA, tibia Imax, tibia Imin, and tibia J (tibia midshaft diameter); QTL on chromosome 2 for femur CA, femur Imax, and femur J (femur length) and for femur Imin (femur midshaft diameter); and QTL on chromosome 10 for tibia CA (tibia length and proximal tibia width).The shared QTL on chromosomes 7 and 2 are also recovered when mapping PC1 in the current study, suggesting they could be determinants of bone morphology across the skeleton.
Several caveats accompany our conclusions.Direct measurement of bone strength ultimately requires mechanical tests (Adams and Ackert-Bicknell 2015), which we did not conduct.Although biomechanical models and empirical comparisons in mice connect the CA and moments of inertia we analyzed with load resistance (Turner et al. 2001;Wergedal et al. 2006;Jepsen 2009;Jepsen et al. 2015), these traits still represent indirect measures of performance.More realistically, cross-sectional measurements capture load resistance abilities rather than performance abilities used by mice as part of their behaviors.Mapping bone strength estimated by mechanical tests could yield different QTL.Nevertheless, our results provide genetic insights into the evolution of determinants of bone strength relevant to load resistance.
Another limitation of our study is its restriction to a single time point, 16 weeks of age.Based on the developmental dynamics of bone structure, we might expect the underlying genetic architecture to depend on age.For example, femur geometry changes with age, and peak bone strength is not reached until 20 weeks in mice from one classical inbred strain (Brodt et al. 1999).
A final caveat is that our evolutionary interpretation is constrained by a lack of knowledge about the most appropriate strain to compare to GI mice.For a mainland reference strain, we used WSB, a wild-derived inbred strain developed by breeding mice trapped in Maryland.Mice are thought to have colonized GI from Western Europe, but the geographic source has proven difficult to resolve (Gray et al. 2014).More generally, some QTL we identified could represent mutations that accumulated along the lineage leading to WSB mice rather than mutations along the lineage leading to GI mice.This inability to polarize QTL is an inherent limitation of all genetic mapping studies that only consider 2 strains.
Our investigation establishes GI mice as a tractable genetic system for understanding how vertebrates evolve to resist loads in their skeletons.Our QTL maps constitute the first step toward identifying specific genes and mutations responsible for the evolution of bone structure in island organisms with unusually large bodies.With the genetic tools available for house mice, the GI mice system is poised to deliver insights into the molecular, cellular, and developmental mechanisms that enable expansions of body size in nature.

Fig. 1 .
Fig. 1.Structural measurements of bone.A micro-CT image from a cross-section of the femur of a GI mouse is shown.Cortical area is estimated as the summed areas of bone in cross-section (lighter shade).Lines indicate orientations for calculating maximum area moment of inertia (Imax) and minimum area moment of inertia (Imin).Polar moment (J) is computed as Imax + Imin.

Table 1 .
Effects of strain and sex on bone structure.

Table 2 .
QTL positions and effects.
expanded bone cross-sectional area and increased moments of inertia without enhancing bone density.This result suggests that the increased skeletal loads imposed by extra weight can be mitigated by structural changes without altering the material properties of bone in an evolutionary context.Structural determinants of load resistance abilities evolved in a coordinated fashion across the bones of GI mice.Correlations

Table 3 .
Candidate genes for QTL that affect the structures of several bones.