Intensive terrestrial or marine locomotor strategies are associated with inter‐ and intra‐limb bone functional adaptation in living female athletes

Abstract Objectives To systematically characterize intra‐limb patterns of skeletal plasticity to loading among living women, in order to better understand regional complexity in structural adaptation within the lower limb and more accurately infer behavior in the past. Materials and methods We used peripheral quantitative computed tomography imaging of the femur, tibia, first and second metatarsals to quantify bone morphology among female controls and athletes representative of either terrestrial or marine mobility, grouped by loading category (odd‐impact, repetitive low‐impact, and high‐magnitude). Parameters included midshaft bone density, areas, rigidity, and shape, epiphyseal bone densities and areas. We assessed between‐group differences and the influence of training history on significant variation among the loading groups. Results Terrestrial mobility strategies were best distinguished by significant midshaft periosteal hypertrophy across the lower limb/foot relative to controls, and by particularly high midshaft femoral and tibial cortical bone areas relative to rowers. Enhanced midshaft bone area was typically paired with decreased bone density among athlete groups. Sport‐specific variation in training duration/timing was significantly correlated with multiple midshaft parameters. Discussion Results demonstrate characteristic patterns of intra‐limb adaptation to terrestrial and marine mobility strategies among active women relative to controls, and highlight components of these patterns that may be shaped in part by differences in loading duration/timing. Additionally, our findings support constraints on skeletal variation in the distal tibia and foot relative to more proximal locations about the knee among living women. For example, metatarsal variation was constrained, but where present reflected sport‐specific variation in force distribution in the foot.

Variation in bone morphology throughout the human lower limb reflects in part competing requirements for structural competency and energetic efficiency. As the cost of accelerating a given limb segment during locomotion is proportional to its mass × moment arm 2 (Hildebrand, 1985), the energetic cost of bone mass is highest most distally in the limb. Theoretically, we might thus predict that bone mass is constrained in distal limb elements relative to more proximal ones, and bone mass distribution in nonhuman mammals does appear to support these predictions. Many mammals, particularly cursorial mammals, exhibit "limb tapering," in which distal bones of the limb are thinner and have lower bone mass than more proximal bones (Currey, 1984;Hildebrand, 1985;Lieberman & Pearson, 2001). Relative to other primates, humans have very long lower limbs for a given body mass, proportions that have significant consequences for our locomotor efficiency. Relatively long lower limbs reduce the metabolic cost of both walking and running (Steudel-Numbers & Tilkens, 2004;Steudel-Numbers, Weaver, & Wall-Scheffler, 2007), but this cost is heavily dependent on the distribution of weight in the limb. For example, moving 3.6 kg of weight from the thigh to the ankle in trained runners increases the metabolic cost of running by 15% (Myers & Steudel, 1985). This is why limb tapering is thought to provide substantial energy savings during the leg swing, as a given amount of mass is less costly to move if it is distributed proximally. This selection for tissue economy has been implicated as a driving force behind variation in bone morphology and robusticity within the limbs among humans (Shaw et al., 2014;Stock, 2006). However, the reduced energy expenditure of limb tapering during locomotion is achieved at the expense of mechanical strength: tapered distal bone has less cortical area through which to dissipate compressive strain, and its reduced periosteal expansion and smaller total area limits its ability to resist bending and torsional strain (Lieberman, Pearson, Polk, Demes, & Crompton, 2003). As a result, tapered distal limb elements experience higher strains than proximal elements, resulting in lower safety factors (strength relative to peak stress during use; Biewener, Thomason, & Lanyon, 1988) and potentially increasing the risk of fracture (Vaughan & Mason, 1975). Quadrupedal mammals, including various primate species, have evolved bone structural adaptations that help reduce the bending moments exerted upon their distal limb segments, by adjusting the orientation of the distal element (Biewener, 1983;Biewener et al., 1988;Polk, 2002), shortening it (Alexander, 1977;Jungers, 1985), and/or reducing its curvature (Bertram & Biewener, 1992;Biewener, 1983;Biewener et al., 1988). The extent to which these relationships between limb tapering, tissue economy, and structural strength extend to the human limb has been considered (Lieberman et al., 2003;Stock, 2006) but not yet explored among living humans. Some cursorial mammals and primate species increase their bone strength relative to mass in distal elements through enhanced intracortical remodeling (Carlson & Patel, 2006;Lieberman et al., 2003;Lieberman & Crompton, 1998;Lieberman & Pearson, 2001;Skedros, Sybrowsky, Parry, & Bloebaum, 2003). This allows strain-induced microdamage to be repaired without increasing bone mass in regions of the skeleton where hypertrophy is energetically costly. Elevated intracortical remodeling among living women has been documented at the tibial midshaft among longdistance runners relative to control subjects (Wilks et al., 2009).
Among archaeological populations of humans, tibial diaphyseal structural properties often appear to reflect inferred locomotion more clearly than femoral diaphyseal properties Macintosh et al., 2014b;Stock, 2006), though see Stock & Pfeiffer, 2001). A similar pattern has been documented among other mammal species: Lieberman and Pearson (2001) found that the tibia exhibited the greatest percent differences in midshaft bending rigidity between exercised sheep relative to controls, compared to both the femur and metatarsal. These patterns appear despite the fact that we might expect the mechanical strain exerted on long bones to increase progressively, proximally to distally, within the limb. In reality, the distribution of strain in the lower limb is regionally variable, affected by body breadth  as well as by characteristics of the gait, limb posture, bone curvature, cross-sectional shape, and soft tissue biomechanics (Skedros et al., 2003). This means that strain attributable to locomotion can be higher in distal limb elements than proximal ones (Stock, 2006), driving structural responses distally despite the associated energetic costs. Whether or not this pattern exists across the femur, tibia, and metatarsals of living humans is unknown, but its characterization would be useful for targeting regions of the lower limb in which adaptation to locomotion may be particularly visible.
The attribution of variation in the most distal limb elements, the hands and the feet, to locomotion is particularly important for understanding human evolutionary history and the transition from arboreality to bipedalism. However, constraints on bone mass and structure should be high here, not just due to the distal position in the limb but also to the functional significance of the hands and feet. Variation among arboreal and terrestrial hominoids in the structural characteristics of cortical bone (metatarsals and/or metacarpals; Jashashvili, Dowdeswell, Lebrun, & Carlson, 2015;Marchi, 2005;Tsegai et al., 2017), trabecular bone (phalanges, metacarpals, and/or metatarsals; Chirchir, Zeininger, et al., 2017;Griffin et al., 2010;Matarazzo, 2015;Stephens et al., 2016;Tsegai et al., 2013), and combinations of the two and bone curvature (Kivell, 2016) in the hands and feet have demonstrated correspondence with locomotor mode. These relationships have been helpful in interpreting hominin and more recent human locomotion in some instances (Dowdeswell et al., 2017). However, the extent to which metatarsal midshaft morphology can inform upon more subtle locomotor differences among habitually bipedal human groups is unclear.
In archaeological human populations, first and/or second metatarsal midshaft CSG properties in terrestrially mobile groups do not vary in response to subtle variation such as terrain complexity (Griffin, Gordon, Richmond, & Anton, 2008), and differ little between prehistoric Jomon hunter-gatherers and modern humans (Hagihara & Nara, 2018). However, male MT1 CSG properties have successfully differentiated between prehistoric human groups with terrestrial relative to marine mobility strategies (Stock & Pfeiffer, 2001), whereas female properties have not. Though this may reflect sex differences in behavior, average size-standardized axial and bending strengths at cadaveric metatarsal midshafts are consistently lower among women than among men (Griffin & Richmond, 2005), and female bone exhibits less responsiveness to mechanical loading than does male bone (Haapasalo et al., 1996;Järvinen et al., 2003;Kannus et al., 1995;Macintosh et al., 2017). Thus, the extent to which diaphyseal adaptation in the metatarsals can inform on even widely varying behavioral strategies among human groups, particularly among women, is uncertain.
The current study uses peripheral quantitative computed tomography (pQCT) scanning to systematically characterize intra-limb patterns of plasticity in response to loading in the femur, tibia, and first and second metatarsals of living women. Athletes from sports involving either terrestrial or marine mobility were grouped into three loading categories: i) odd-impact loading (soccer), ii) repetitive low-impact loading (running), and iii) high-magnitude loading (high load magnitudes exerted by muscle, but no ground impact; rowing). These loading groups are described in detail below. Between-group differences relative to controls were assessed in midshaft bone density, area, rigidity, and shape, as well as epiphyseal bone densities and areas. We also assessed the influence of both the duration of sport participation and its timing relative to menarche on any significant variation in bone outcomes among the sport groups relative to each other. By characterizing patterns of cortical and trabecular bone functional adaptation in the lower limb among living women, we aimed to test the correspondence between loading and the properties and regions of the limb typically utilized in anthropological research to interpret behavior. In doing so, we hope to enable the more nuanced and confident interpretation of behavior in the past, particularly among women.

| Participants
Sports were selected for inclusion in the study that reflected both intensive terrestrial and marine mobility strategies, subdivided into loading categories defined by Nikander et al. (2006, and Nikander, Sievänen, Heinonen, and Kannus (2005).
Intensive terrestrial mobility strategies were represented by two sport groups, both of which involve loading exerted through ground impacts: odd-impact loading (soccer; N = 11) and repetitive lowimpact loading (endurance running; N = 17). Odd-impact loading involves ground impacts associated with rapid turning, acceleration and deceleration, and a mix of sprinting and running. The directionality of this loading was considered atypical; its rapid changes of direction, engendering large torque forces, would not typically be encountered in everyday terrestrial locomotion (Niinimäki et al., 2017). Our other terrestrial loading group, repetitive low-impact loading, involves lowintensity ground impacts applied with relatively constant speed over a long period of time. The directionality of this loading was considered typical of everyday terrestrial locomotion, being predominantly in the antero-posterior direction and not involving rapid changes of direction or high torque. Marine mobility strategies were represented by one loading group, in which ground impact is absent from sport-specific movements: high-magnitude loading (sweep rowing; N = 17). Highmagnitude loading does not exert ground impact loads, but rather exerts high load magnitudes through maximal muscle forces, applied through coordinated movements. In this case, among rowers, the loading involves no rapid changes of direction. In sweep-style rowing, the athletes seat moves while their feet are attached to a static footplate. With each stroke, the athlete drives off of this footplate with powerful muscular contractions of the lower limbs, producing high muscle magnitudes and muscle joint contact forces acting on the foot, ankle, and knee (Hase et al., 2002), but no vertical ground reaction forces or ground impacts are experienced by the lower limbs. Athletes were compared to a group of recreationally active control subjects (N = 26).
All participants were healthy premenopausal adults of European descent living in the United Kingdom, and were aged 19-43 years.
We did not control for handedness in our study, as all participants were right-handed with the exception of three control subjects and one rower. The following exclusion criteria were applied to all subjects: any medical condition or medication known to interfere with bone metabolism, any current or recent (past 12 months) pregnancy or lactation, 18 years of age or younger, or postmenopausal status.
Additional exclusion criteria for athletes were: participation in the sport of interest for fewer than 3 years, any significant injury within the past year that rendered them inactive for over 1 month, and any current intensive participation in another sport other than the one for which they were recruited. Additional exclusion criteria for control subjects were: any current or past participation in competitive sport and any current or past participation of more than 3 hr a week of weight-bearing intensive physical activity.

| Anthropometrics
Height was measured to the nearest 0.01 cm using a SECA 274 stadiometer. Weight was recorded to the nearest 0.1 kg with a SECA electronic scale. Femoral length and maximum tibial and first metatarsal (MT1) lengths were obtained from participants using sliding calipers according to the methods in International Standards for Anthropometric Assessment (International Society for the Advancement of Kinanthropometry, 2001). Femoral length was measured as the distance between the proximal border of the greater trochanter and the distal border of the lateral condyle, so does not represent the true maximum length of the bone. Tibial length was measured as the distance between the proximal medial border of the tibial plateau and distal border of the medial condyle. First metatarsal length was measured as the distance between the proximal and distal joint lines.

| Peripheral quantitative computed tomography and properties of interest
All cross-sectional bone images were collected using peripheral quantitative computed tomography (XCT-3000; Stratec Medizintechnik GmbH, Pforzheim, Germany) at the PAVE Imaging and Performance Laboratory in the Department of Archaeology at the University of Cambridge. Cross-sectional images were obtained at 50% and 4% of maximum length (from the distal end) in both the right femur and tibia, and at 50% of length in the right first metatarsal (MT1). All second metatarsal (MT2) properties were quantified from the section location visible in the pQCT image taken at the MT1 midshaft, a location considered to approximate the MT2 midshaft. Because it was not possible to obtain true maximum length from the femur, the``midshaft" slice used here is slightly distal to true midshaft, though is still taken close to diaphyseal transverse cross-sectional minima. Movement artifacts affecting bone contours were present in five femoral midshaft slices, three midshaft MT1 slices, and four midshaft MT2 slices. Data from these affected slices were removed prior to data analyses, and sample sizes per group and per bone are provided in the legends for Tables 2 and S1. To assess cross-sectional bone strength and shape, crosssectional images of midshaft femora, tibiae, and metatarsals were imported into ImageJ (http://rsbweb.nih.gov/ij/) and parameters were quantified using BoneJ, a bone image analysis plug-in (Doube et al., 2010). For femora and tibiae, the "Optimize Threshold" function was used to isolate cortical bone. For metatarsals, pixels with 8-bit brightness between 128 and 255 were considered cortical bone. To assess bending/torsional rigidity, polar second moments of area (J; mm 4 ) were calculated from the midshaft femur, tibia, MT1, and MT2. The polar second moment of area provides a measure of twice average bending and torsional rigidity of the cross-section, and is the sum of any two perpendicular second moments of area (Ruff, 2008); in this case, second moments of area about the maximum (I max ; mm 4 ) and minimum (I min ; mm 4 ) axes of the cross-section were used. To assess cross-sectional shape, a ratio of I max /I min was calculated, which provides information on the distribution of bone in cross-section along these maximum and minimum axes.

| Statistical analysis
Data were assessed for normality using the Kolmogorov-Smirnov test and assessments of skew and standard error. Variables that were nonnormally distributed were natural log transformed prior to analyses. Group differences in anthropometric variables, and athlete differences in training variables, were assessed using one-way analysis of variance (ANOVA), using Hochberg's GT2 or Games-Howell post-hoc tests. Between-group differences in bone outcomes were assessed by analysis of covariance (ANCOVA) using age, weight, and height as covariates. No confidence interval adjustment was used due to small sample sizes. Percentage differences for adjusted means and 95% CIs between athlete groups relative to the reference group (control subjects) were calculated using age-, height-, and weight-adjusted means from ANCOVA, and antilogged where required. Bivariate Pearson's correlations were used to assess the relationships between training history variables, age, height, weight and the bone parameters (unadjusted) that differed significantly among the athlete groups in ANCOVA. If age, height, and/or weight were significantly correlated with any of these bone parameters, these were subsequently controlled for in partial correlations, to determine if the relationship between training history and bone parameters remained. All analyses were performed using SPSS version 23.0. Alpha was set as p < .05 for all analyses.

| RESULTS
Descriptive statistics and group differences for anthropometric and training history variables (ANOVA) are given in Table 1. Comparisons of bone outcomes adjusted for age, height, and weight are given in Table 2 (ANCOVA), with unadjusted and raw antilogged summary statistics for each group provided in the Supplementary Information (Table S1). Percent differences and 95% CIs for adjusted marginal means by bone property among athlete groups relative to the reference group (controls) are presented in Figure 1

| Between-group differences: femur
At the femoral midshaft, the odd-impact loading group (soccer) exhibited the most widespread adaptation relative to controls. Soccer players had significantly greater midshaft TA, J and I max /I min (by~9 to 19%) with no significant change in MA, suggestive of hypertrophy at the periosteal surface. Soccer players also had significantly higher CA (~13%) paired with significantly lower CBD (~2%) than control subjects. The high-magnitude loading group (rowers), lacking ground impact, did not exhibit any significant change from controls in variables reflecting periosteal hypertrophy. However, they did exhibit a pattern of significantly greater femoral CA (~6%) paired with significantly lower CBD (~2%) than controls.
The repetitive low-impact group (runners) significantly exceeded control subjects in midshaft femoral TA and J (by~10 to 20%) but not MA, suggesting periosteal hypertrophy. Runners also significantly exceeded controls in midshaft CA (~14%), but this was not paired with any significant corresponding change in CBD.
At the distal femoral epiphysis, adaptation among athletes was consistent across all loading groups relative to controls. In the cortical bone shell, all exercise loading groups exhibited significantly higher adjusted mean CA (by~14 to 18%) and CBD (by~10 to 14%) than controls. In the trabecular bone compartment, all loading groups exhibited significantly greater adjusted mean TrabBD than controls (by~8 to 13%). There were no significant differences in TrabA at the distal femur between any of the groups.
Among the athletes, the groups performing intensive terrestrial locomotion (running and soccer) differed significantly from each other in midshaft femoral CBD, but both groups exhibited significantly higher midshaft femoral CA than our marine mobility group, rowers. Soccer players also exceeded rowers in midshaft femoral J. There were no significant differences among athlete groups at the distal femoral epiphysis.

| Between-group differences: tibia
At the tibial midshaft, the loading groups involving intensive terrestrial locomotion (soccer and endurance running) exhibited significantly enhanced TA and CA (by~9 to 15%) relative to control subjects, with no significant change in MA. The repetitive, low-impact loading of endurance running was also associated with significantly enhanced J and I max /I min (by~20 to 23%) relative to controls. The odd-impact loading of soccer was associated with significantly enhanced J as well, by~18% relative to controls, but no significant change in I max /I min . All athlete groups exhibited lower FIGURE 1 Grouped scatterplot for percent differences and 95% confidence intervals in adjusted marginal means of bone outcomes among athletes relative to controls at the distal femoral and tibial epiphyses (4% section location). Fem = femur; Tib = tibia. The 0% vertical dotted line indicates the mean of the control subjects. * indicates significantly different than control subjects FIGURE 2 Grouped scatterplot for percent differences and 95% confidence intervals in adjusted marginal means of bone outcomes among athletes relative to controls at the femoral and tibial midshaft (50% section location). Fem = femur; Tib = tibia. The 0% vertical dotted line indicates the mean of the control subjects. * indicates significantly different than control subjects adjusted mean tibial midshaft CBD than controls, but these differences were only significant among soccer players and rowers (~1 to 2%).
At the distal tibial epiphysis, all athlete groups exhibited significantly higher adjusted mean CA (by~10 to 15%) and mean TrabBD (by~10-14%) than controls. Only soccer players exhibited lower distal CBD than controls (by~7%), while distal CBD among rowers significantly exceeded controls (by~9%). The repetitive, low-impact loading of endurance runners was instead associated with enhanced TrabA relative to controls (by~8%).
Among the athlete groups themselves, the tibiae of groups performing intensive terrestrial locomotion (soccer players and runners) differed significantly from tibiae of the marine mobility group (rowers) in several instances. Soccer players had significantly higher midshaft tibial CA than rowers, while endurance runners had significantly higher midshaft tibial J and distal tibial TrabA than rowers.
Further, endurance running was associated with significantly higher tibial I max /I min (more elliptical in cross-section) than both soccer or rowing.

| Between-group differences: metatarsals
At the metatarsal midshafts, significant between-group differences were the fewest of all lower limb section locations, and were only documented among the groups involving intensive terrestrial locomotion. Endurance runners had significantly lower MT1 I max /I min (more circular in cross-section) than controls (by~13%), while soccer players had significantly higher MT2 J than controls (by over 28%). Among athlete groups, the low I max /I min of the MT1 midshaft of runners was significantly exceeded by both soccer players and rowers.

| Impact of training history on bone outcomes among athletes
All training history variables differed significantly among athlete groups (see Table 1). Soccer players began training significantly earlier relative to menarche than runners and rowers, and had trained for significantly more years than rowers. Despite this, their current training volume was relatively low, and was significantly exceeded in hours per week by both rowers and runners. Training volume was particularly high in the rowing group, significantly so relative to both other athlete groups.

| DISCUSSION
The characterization of specific patterns in cortical and trabecular bone functional adaptation throughout the lower limb and foot of living humans with known loading is crucial in order to more accurately identify skeletal correlates of behavior in the past. A summary of the main findings of the current study, by loading group, is presented in Table 4. By assessing bone mass hypertrophy and structural variation in response to three different terrestrial and marine loading regimes (repetitive low-impact loading, odd-impact loading, high magnitude loading) among living women, we document evidence of substantial midshaft adaptation to loading in the tibia, followed closely by the femur, but markedly less adaptation in the metatarsal midshafts. This intra-limb pattern of midshaft adaptation, whereby the tibia can demonstrate stronger responsiveness to loading than the femur despite its more distal location in the limb, supports patterns documented among human skeletal remains with inferred behavior Shaw et al., 2014) and from other mammals (Lieberman & Pearson, 2001), likely relating to regional variability in the distribution of strain throughout the limb in locomotion. Further, results demonstrate that bone mass hypertrophy and structural variation in the most distal regions of the limb, the ankle and foot, is constrained relative to more proximal locations about the knee (midfemur through midtibia) in living women. This finding supports theoretical predictions based on trade-offs between bone strength and energetic efficiency acting on these regions in particular. Where significant sport-specific patterns of distal limb adaptation were documented, we highlight their functional relevance to the interpretation of terrestrial locomotion from these constrained regions.

| Lower limb bone adaptation to terrestrial and marine mobility strategies
The loading experienced by a limb bone during terrestrial locomotion varies within the limb but also within a bone itself. For example, within the tibia during walking, loading at the distal end is primarily axial, whereas as section location moves proximally, axial strain interacts with the tibias diaphyseal curvature to create bending and torsional loads that increase progressively with increasing proximity to the knee (Bertram & Biewener, 1988;Biewener, 1991;Capozza et al., 2010;Garcia & da Silva, 2004;Wehner, Claes, & Simon, 2009  locomotion (soccer and running) that differed significantly from controls in midshaft tibial TA, J, and I max /I min . These midshaft adaptations were accompanied by simultaneous distal adaptations in the tibial epiphysis that increase compressive strength while minimizing mass, via enhanced TrabBD. Due to the nonlinear relationship between density and strength in trabecular bone (Ebbesen et al., 1997;Rittweger et al., 2006), small increases in density produce large increases in compressive strength (Wilks et al., 2009). The odd-impact loading of soccer, with its rapid changes of direction and speed, was associated with particularly high TrabBD in the distal tibia, exceeding controls by 14%.
In contrast, the marine mobility group, comprised of rowers, never differed significantly from controls in midshaft variables reflecting this diaphyseal adaptation to the bending/torsional loads (TA, J, I max /I min ) characteristic of intensive terrestrial locomotion. However, rowers did exhibit significant adaptation to axial strain relative to controls, via enhanced midshaft CA and epiphyseal TrabBD and CA relative to controls in both the femur and tibia. This may be attributable to the joint contact forces exerted on the lower limbs during sweep rowing. During the recovery phase of the rowing stroke, the athlete compresses the lower limb as they move forward toward the foot plate, flexing the ankle joint and eliciting up to 2,800 N of compressive joint contact force (Hase et al., 2002). The athlete then drives off the footplate with powerful muscular contractions that extend the legs, causing up to 370 N of vertical force on the foot against the plate and exerting peak joint contact forces in compression of 4,100 N on the knee (Hase et al., 2002); these forces exceed 6× the rowers body weight and are twice as high as those experienced during walking or cycling (2-3× body weight; Anderson & Pandy, 2001;Neptune & Kautz, 2000). The marine mobility strategies employed by prehistoric populations would have involved boats with fixed seats, like canoes or kayaks, so individuals would not have been driving off of a footplate.
Future research involving limb bone adaptation among modern canoe or kayak athletes would be beneficial to provide a more accurate assessment of the effects of traditional marine mobility on the skeleton.

| Directionality among terrestrial loading regimes
Midshaft lower limb bone shape ratios differed significantly between our intensive terrestrial loading groups in two instances. The directionality of loading in endurance running, predominantly anteroposterior and lacking rapid turning movements, was associated with significantly higher midshaft tibial I max /I min and lower MT1 I max /I min than soccer players. An extreme example of the pattern of tibial change characteristic of endurance runners is presented in Figure 4.
The enhanced I max /I min in the shank is likely reflecting relative anteroposterior hypertrophy at the periosteal surface among endurance runners, as evidenced by the tibia in Figure 4a. Similarly, the lower mean I max /I min at MT1 likely reflects enhanced dorsoplantar diameters.
However, runners did not differ significantly from soccer players in femoral I max /I min , as the latter exhibited particularly high values here, sufficient to significantly exceed control subjects (see Figure 5). The rapid changing of direction and speed in soccer results in large torques about the hip, as the gluteal muscles contract to extend the hip as the athlete moves and turns (Niinimäki et al., 2017). It is possible that the muscle activity required to stabilize the hip and provides powerful hip extension during acceleration could be contributing to the shape change among soccer players in the femur relative to more distal locations in their limb, and relative to control subjects. However, Niinimäki et al. (2017) did not find any significant changes in shape ratios at the proximal or midshaft femur among soccer players relative to other sporting groups or to controls in their study of female athletes.

| Evidence of midshaft intracortical remodeling in the lower limb in response to strain
The significantly enhanced midshaft femoral and tibial CA among all athlete groups relative to controls was consistently paired with significantly lower CBD at these locations among both soccer players and rowers (but not runners). Lower CBD at the tibial midshaft among athletes relative to controls has been reported previously in the literature (Liu et al., 2003;Rantalainen et al., 2011;Wilks et al., 2009), and is often highly correlated with increased porosity (Bousson et al., 2000). Thus, reduced midshaft CBD relative to controls in some loading groups may reflect enhanced intracortical remodeling rates in response to microdamage from repeated loading cycles (Burr, Martin, Schaffler, & Radin, 1985;Mori & Burr, 1993). If so, our results provide evidence for intracortical remodeling in response to mechanical strain not just at the tibial midshaft among living women, but at the femoral midshaft as well. It is unclear why runners do not follow this pattern.
However, midshaft CBD exhibited a significant positive relationship with the timing of training relative to menarche (see below); lower CBD was associated with earlier initiation of sport relative to menarche, and endurance runners were the group that initiated their training the latest relative to menarche (see Table 1).  (Doershuk et al., 2019;Mittra et al., 2005;Saers et al., 2016). For example, mean trabecular thickness in the calcanei of male endurance runners is significantly positively related to their weekly running distance and years of running (Best et al., 2017). The comparison of two-and threedimensional trabecular bone imaging in the future may enable the more nuanced interpretation of variation in pQCT-derived trabecular bone variables among living women. In any case, the maintenance of significant bone mass hypertrophy so distally in the limb solely among endurance runners suggests that this bone mass is necessary to accommodate the specific strains associated with the low-impact repetitive loading of running, despite the energetic costs.

| The functional significance of distal limb bone structural variation: metatarsal midshafts
The significantly more "circular" (low I max /I min ) MT1 midshafts of endurance running group relative to all others likely reflects enhanced dorsoplantar relative to mediolateral dimensions; the former were typically smaller than the latter among the women in this study (data not shown). This relative dorsoplantar expansion at MT1 may be related to the particularly high peak pressures that running exerts on the MT1 region, which increase by up to 72% relative to those exerted during walking (Rozema, Ulbrecht, Rammer, & Cavanagh, 1996). Similarly, the flexor tendons of the longitudinal arch load MT1 heavily during locomotion, resulting in forces on the first metatarsal head of up to 119% of body weight (Jacob, 2001 The atypical turning movements of soccer, such as "cutting", exert much higher plantar pressures on the foot than does running in a straight line (Eils et al., 2004;Orendurff et al., 2008), engendering substantial bending moments on lateral regions of the foot that are not typically accustomed to these (Eils et al., 2004;Ekstrand & Torstveit, 2012). This particular pattern of metatarsal loading may explain why the odd-impact loading of soccer players was associated with uniquely enlarged and strengthened MT2 midshafts both relative to the other loading groups and to their MT1, as seen in Figure 6. In addition, soccer players often kick the ball with the dorsolateral aspect of the foot, a further source of strain on MT2. This metatarsal exhibits the weakest CSG properties relative to its length of all metatarsals (Griffin & Richmond, 2005) and it is the only metatarsal that does not demonstrate a large safety factor between bone structure and peak force (Lidtke, Patel, & Muehleman, 2000). Thus, the significantly enlarged MT2 J among soccer players relative to controls (28% higher) may reflect adaptation to reduce this mismatch between midshaft MT2 cross-sectional robusticity and the particularly high mechanical loads exerted upon this region in soccer.
4.5 | Training history affects between-group differences in some lower limb variables In our study, the timing of training relative to menarche was significantly positively associated with femoral CBD (r = 0.383, p < .05): lower CBD was associated with earlier initiation of training relative to menarche. Thus, the significantly lower femoral midshaft CBD of soccer players relative to the endurance running group may be reflecting, in part, the prepubertal sport specialization of the former (see Table 1). However, CBD is not typically assessed in anthropological studies of prehistoric behavior, while femoral and tibial bone areas and bending/torsional rigidity are; these were occasionally influenced by the overall duration of training.
In our study, the total number of sport-specific years of training was significantly associated with femoral and tibial midshaft J and tibial midshaft CA (r = 0.529-0.615, p < .001). Despite differing significantly in the timing of their training relative to menarche, soccer players and endurance runners did not differ significantly in the total duration of that training (between 10 and 12 years), and they did not differ in mean CA or J. In contrast, soccer players did differ significantly from the rowing group in the total duration of training (having trained an average of 5.5 years longer); these two groups differed significantly in midshaft femoral J and tibial CA. Thus, differences in midshaft femoral and tibial FIGURE 5 Right femoral midshaft pQCT images demonstrating cortical bone shape adaptation across loading groups. Each athlete exhibits I max /I min closest to the mean for their group. Relative to the control subject, the soccer player exhibits 12% higher I max /I min , the runner exhibits 8% higher I max /I min , and the rower exhibits 3.5% higher I max /I min . Images not to scale. Anthropometric characteristics as follows: Soccer player-age: 26 yrs; height: 167.1 cms; weight: 65.1 kgs. Endurance runner-age: 30 yrs; height: 160.0 cms; weight: 50.7 kgs. Rower-age: 22 yrs; height: 184.7 cms; weight: 75.7 kgs. Control subject-age: 20 yrs; height: 159.0 cms; weight: 57.1 kgs. Training characteristics as follows: Soccer player-age at menarche: 14; sport-specific years: 18; current hrs/wk: 6. Endurance runner-age at menarche: 10; sportspecific years: 7; current hrs/wk: 10.0. Rower-age at menarche: 12; sport-specific years: 5; current hrs/wk: 20.0 FIGURE 6 Midshaft pQCT images of right MT1 and MT2 by loading group demonstrating percent differences in bending/torsional rigidity (J) between the first and second metatarsals CA and J between loading groups may be reflecting not just the type of loading, but the cumulative duration of loading as well.
The fact that the timing of loading relative to puberty was not significantly correlated with midshaft femoral and tibial CA in particular is unexpected, as variation in timing has been shown to have strong effects on percent differences in TA and CA between the upper limbs among females in racquet sports (Bass et al., 2002;Haapasalo et al., 1996;Kontulainen et al., 2002). For example, Haapasalo et al. (1996) showed significantly larger side-to-side differences in humeral cortical wall thickness among young starters (began at~9 years of age) than old starters (began at~29 years of age) despite twice as many training years in the latter (mean of 8 and 14 total years of training, respectively). However, the interaction of training duration and its timing is complex; among male racquet-sport athletes, accounting for training duration eliminated significant differences in humeral TA between pre-and peri-pubertal players (Ducher, Daly, & Bass, 2009). Similarly, in the tibia, elite young adult male soccer players who had been training for more than 3 years had significantly higher TA and CA than players with fewer years of training (Hart et al., 2016). Our results suggest that more work is needed to understand the interaction between loading duration and timing relative to menarche on bone structure among women in particular.
In prehistoric and archaeological populations, it is often very difficult to account for the influence of time/timing on bone functional adaptation, as we simply do not know the age at which individuals initiated adult behaviors, particularly relative to puberty, or for how long they may have been participating in these behaviors. If the cumulative duration of loading is indeed contributing to variation in CSG properties such as CA or J in the lower limb, then this has implications for the interpretation of behavior from these properties when comparing groups with, for example, very different age profiles or different degrees of individual task specialization. Though the patterning of loading bouts throughout the day and the lifetime may differ between prehistoric women and modern-day athletes, the identification of relationships between patterns of bone structural variation and known loading characteristics and history is a significant step toward the more accurate interpretation of skeletal variation in the past among women.

| CONCLUSION
By systematically assessing between-group differences in femoral, tibial, and first and second metatarsal bone parameters, we provide evidence of an optimized balance between bone mass and safety factors in the lower limb of living women, with greater constraints on bone mass hypertrophy and structural variation in the distal tibia and foot relative to more proximal locations about the knee (midfemur through midtibia). With regards to the interpretation of mobility strategy among past human populations, our results highlight characteristic patterns of intra-limb adaptation to intensive terrestrial mobility strategies among active women relative to controls and to athletes employing a marine mobility strategy. By characterizing these intralimb patterns of loading-related structural variation among living women, we indicate the combinations of properties and regions of the limb that best reflect specific loading characteristics among living women, and highlight the importance of considering intra-limb location and loading history when attempting to interpret behavior in the past.

ACKNOWLEDGMENTS
The authors would like to thank Dr. Peter Ellison and the two reviewers whose comments were helpful in improving this manuscript. We would also like to thank all of the women who volunteered