Computed tomography shows high fracture prevalence among physically active forager-horticulturalists with high fertility

Modern humans have more fragile skeletons than other hominins, which may result from physical inactivity. Here, we test whether reproductive effort also compromises bone strength, by measuring using computed tomography thoracic vertebral bone mineral density (BMD) and fracture prevalence among physically active Tsimane forager-horticulturalists. Earlier onset of reproduction and shorter interbirth intervals are associated with reduced BMD for women. Tsimane BMD is lower versus Americans, but only for women, contrary to simple predictions relying on inactivity to explain skeletal fragility. Minimal BMD differences exist between Tsimane and American men, suggesting that systemic factors other than fertility (e.g. diet) do not easily explain Tsimane women’s lower BMD. Tsimane fracture prevalence is also higher versus Americans. Lower BMD increases Tsimane fracture risk, but only for women, suggesting a role of weak bone in women’s fracture etiology. Our results highlight the role of sex-specific mechanisms underlying skeletal fragility that operate long before menopause.


Introduction
Compared to other hominoids and extinct hominins, modern humans have postcranial skeletons that are more gracile (i.e. lower bone mass and strength for body size) (Cotter et al., 2011;Ruff et al., 2015;Ryan and Shaw, 2015). Declining skeletal strength has been documented in the cortical structure of long bone diaphyses (e.g. size or shape in a cross-section) and in trabecular bone micro-structure (e.g. thickness, bone volume fraction) and is particularly evident in the later Pleistocene or Holocene (Chirchir, 2019;Chirchir et al., 2015;Ruff et al., 2015;Ryan and Shaw, 2015). Physical inactivity is the most common explanation for skeletal fragility, based largely on the well-established principle that impact forces from load bearing and muscle contraction trigger bone deposition (Carter and Orr, 1992). According to this principle, which has been documented in numerous studies of competitive athletes and exercise interventions (e.g. Polidoulis et al., 2012;Warden et al., 2014), and inferred from skeletal remains of prehistoric populations (e.g. Macintosh et al., 2017;Ruff et al., 2015), bone responds to physical activity demands by adding tissue and altering crosssectional distribution in the direction of highest bending strains (i.e. change in length per unit length) (but see Demes et al., 2001;Lieberman et al., 2004;Lovejoy et al., 2003 and references therein). This mechanical response of bone to loading is variable throughout the body, depending on types of weight-bearing activity and muscle function.
Evolutionary life history theory provides a broad explanatory framework that incorporates ultimate and proximate levels of analysis for understanding variability in bone strength (i.e. ability to withstand an applied load). In all organisms, limited resources are allocated to competing metabolic demands so as to optimize biological fitness. Due to higher fitness gains of reproduction earlier versus later in life (Williams, 1957), natural selection often prioritizes investments in earlier reproduction over somatic maintenance. Organisms may thus increase fertility at the expense of maintenance (Kirkwood and Austad, 2000), and we should expect greater energetic investments in reproduction to trade-off against investments in maintenance (Stearns, 1992). Consistent with life history theory, reproductive effort is expected to moderate effects of physical activity on adult bone strength, which is an indicator of energetic investment in somatic maintenance.
Tests of this life history trade-off in humans are inconclusive (Le Bourg, 2007), in part because many studies focus on mortality rather than investments in maintenance per se (but see Ryan et al., 2018;Ziomkiewicz et al., 2016), precluding direct analysis of whether greater reproductive effort inhibits maintenance. Bone tissue is ideal for examining metabolic trade-offs between reproduction and maintenance. Constant remodeling is necessary to maintain bone strength, and the skeleton is a general mineral reservoir for the competing metabolic demands of maternal maintenance and fetal bone accretion or lactation (Stieglitz et al., 2015). The average full-term human fetus has~30 g calcium, 20 g phosphorus and 0.8 g magnesium, and at least 80% of these macro-minerals are accreted in the third trimester (see Kovacs, 2016 for a comprehensive review of bone metabolism during pregnancy, lactation and post-weaning). For an average-sized fetus this corresponds at week 24 of gestation to a mean calcium [phosphorus] transfer rate of~60 mg/day [~40 mg/day] and between weeks 35-40 of 300-350 mg/day [200 mg/day]. In the third trimester, hourly fetal transfers of calcium and phosphorus are between 5 and 10% of that present in maternal plasma, which is enough to provoke maternal hypocalcemia and hypophosphatemia. Generally, patterns of bone turnover are similar comparing pre-to early pregnancy states, but turnover increases during the third trimester to create a net resorptive state. During the first six months [second six months] postpartum,~200 mg/ day [~120 mg/day] of calcium is secreted into human breast milk. Analyses of bone turnover markers, bone mineral density (BMD), and bone structure by high-resolution peripheral quantitative computed tomography (HR-pQCT) suggest that lactating women are in negative calcium balance, especially when milk production is elevated. Longitudinal studies show consistent declines in lactating women's BMD or bone mineral content, with mean declines of 3-10% after 3-6 months of lactation. The greatest BMD losses (5-10%) occur in the lumbar spine, with modest losses (<5%) occurring at sites with less trabecular bone, and smaller losses (<2%) at sites containing mostly cortical bone. The few studies utilizing HR-pQCT in the limbs (radius, femur) show smaller (<2%) reductions in trabecular thickness and cortical thickness and volume, with greater reductions among women who lactate longer. Thus, because mineral allocations to maternal maintenance and reproduction draw from the same skeletal reservoir, direct trade-offs between these competing demands should manifest in bone. These trade-offs are expected to manifest in the longer term, regardless of whether maternal bone tissue fully or only partially recovers from mineral losses following a specific bout of gestation and lactation. Life history theory makes no assumptions or predictions about the extent of bone mineral recovery (i.e. whether full or partial) after weaning a specific child. Accordingly, a general hypothesis from life history theory is that greater reproductive effort constrains the ability of bone tissue to respond to mechanical loading and high physical activity levels (PALs). This hypothesis is not an alternative to and may complement other hypotheses of bone structural variation emphasizing developmental factors affecting the trade-off between investment in growth and reproduction (Macintosh et al., 2018). But unlike other hypotheses of bone structural variation derived from life history theory or proximate explanations (focusing, for example on nutrition, inflammation, hormones), the hypothesis emphasizing effects of reproductive effort uniquely predicts sex differences within and between populations, the magnitude of which should be influenced by relative investment in reproduction.
Timing of reproduction, in addition to lifetime reproductive effort, is also expected to affect bone strength (Macintosh et al., 2018). Peak bone mass is typically not achieved until the late 20s, so earlier pregnancy and lactation can disrupt maternal bone growth and/or mineralization (Madimenos et al., 2012;Stieglitz et al., 2015), potentially reducing peak bone mass and thus later-life bone strength (but see Chantry et al., 2004). In addition to early onset of reproduction, short interbirth intervals (IBIs) can potentially generate unbalanced cycles of maternal bone resorption and formation, limiting maternal skeletal recovery before subsequent pregnancy (Stieglitz et al., 2015). Whether the maternal skeleton is fully or partially restored post-weaning to its prior mineral content and strength, despite potentially lasting micro-architectural changes shown in recent imaging studies, is currently debated (Clarke and Khosla, 2010;Kovacs, 2016;Kovacs, 2017;Wysolmerski, 2010) and physiological mechanisms underlying post-weaning maternal skeletal recovery are not well understood. Dual-energy X-ray absorptiometry data suggest that for most women, lactation-associated BMD losses and micro-architectural deficits are reversed by 12 months post-weaning, although patterns can vary by skeletal site and most studies are conducted among well-nourished women with low fertility (Kovacs, 2016). Restorative capacity partly depends on lactation duration; HR-pQCT studies of the radius and femur show recovery of trabecular and cortical micro-architecture in women lactating for shorter periods, but incomplete recovery in women lactating for longer periods. A recent HR-pQCT study found incomplete recovery of trabecular and cortical micro-architecture in the tibia and radius after a median of 2.6 years post-weaning for women exclusively breastfeeding for five months (Bjørnerem et al., 2017). Incomplete restoration of lumbar spine BMD to pre-pregnancy values by 12 months post-partum is observed among rural Gambian women practicing on-demand breastfeeding for about two years (Jarjou et al., 2010).
A relevant literature on 'maternal depletion syndrome' (Jelliffe and Maddocks, 1964) examines trade-offs between reproductive effort and maternal health more broadly in high fertility contexts. Currently, evidence for maternal depletion is mixed Tracer, 2002), though most prior studies focus only on the period covering one or two births rather than the cumulative long-term effects on health of repeated pregnancies. Moreover, most of these prior studies focus on maternal anthropometric status (e.g. weight, adiposity) rather than bone tissue per se. Bone tissue fluctuates less than anthropometric markers with short-term changes in energy balance, rendering studies of bone less susceptible to sampling biases.
Here, we examine in a natural fertility population, Tsimane forager-horticulturalists of Bolivia, whether greater reproductive effort compromises bone strength, particularly for women given their greater energetic costs of reproduction. Tsimane are an ideal population to test whether women's greater reproductive effort compromises bone strength. Tsimane fertility is high (total fertility rate = 9 births per woman), birth spacing is short (Stieglitz et al., 2015), breastfeeding is ondemand, effective birth control is rare, and PALs are high (Gurven et al., 2013), as is typical of other small-scale rural subsistence populations. In a population-representative sample of adults aged 40+ years, who mostly have completed their reproduction, we utilize thoracic computed tomography (CT) to measure two primary indicators of bone strength in thoracic vertebrae: BMD, which accounts for~70% of the variance in bone strength (NIH Consensus Development Panel on Osteoporosis Prevention, Diagnosis, and Therapy, 2001), and fracture prevalence and severity. We focus on thoracic vertebrae since spontaneous thoracic vertebral fractures are among the most common osteoporosis-related fractures in humans (Sambrook and Cooper, 2006). Such fractures have not been observed in wild or captive apes, even in individuals with severe osteopenia (Gunji et al., 2003), suggesting that modern humans are more susceptible than other primates to osteoporosis-related fractures (Cotter et al., 2011).
Most activities of daily living, including sitting, walking, running and lifting, generate loads on human vertebrae (Myers and Wilson, 1997;Rohlmannt et al., 2001;Stewart and Hall, 2006), and thus thoracic vertebrae track mechanics of both lower and upper limbs. Even regular breathing appears to generate intradiscal pressure and some loading in the thoracic spine (Polga et al., 2004). Vertebral bodies, which are composed mostly of trabecular bone surrounded by a thin cortical shell, function largely as shock absorbers and can deform to a greater degree than tubular bones (Seeman, 2008); deformation facilitates spinal flexion, extension and rotation. The major loading mode on human vertebral bodies is axially compressive, and most axial force is carried by the trabecular bone (Myers and Wilson, 1997). For many activities (e.g. neutral standing, standing with weight, mild trunk flexion and extension, lifting objects above the head), the greatest compressive loads along the spine are generated in the thoracolumbar region (Bruno et al., 2017). The greatest compressive vertebral loads occur during activities in which body mass or externally applied weights are shifted anteriorly, such as during trunk flexion or carrying weight in front of the body (Bruno et al., 2017;Polga et al., 2004;Rohlmannt et al., 2001). To accommodate forces, the architecture of the vertebral body trabecular bone consists of thick vertical plates and columns supported by thinner horizontal trabeculae. This trabecular structure changes with age, such that vertical plates are successively perforated during remodeling and converted into columns, whereas horizontal trabeculae perforate and disappear (Mosekilde et al., 1987). These age-related trabecular structural changes can result in vertebral strength declines that are greater than predicted from bone mass estimation alone. Vertebral strength is compromised more by loss of trabecular connectivity than by trabecular thinning, and women are more susceptible than men to age-related horizontal trabecular perforation and disappearance (Mosekilde, 1989).
In this paper, we first test, among Tsimane women, whether greater reproductive effort -indicated by earlier age at first birth, higher parity and shorter IBIs -is associated with reduced thoracic vertebral BMD. We then test whether Tsimane BMD is lower, particularly for women, than a matched American sample with directly comparable CT-derived indicators of bone strength. This latter prediction follows from the hypothesis that greater reproductive effort compromises bone strength, which is consistent with a trade-off between investment in reproduction and maintenance as posited by evolutionary life history theory. In contrast, a simple prediction from a physical inactivity hypothesis for compromised bone strength posits the opposite, that is, lower BMD for Americans than Tsimane, for both sexes, given lower PALs, on average, among Americans. Regarding the second measured bone strength indicator, thoracic vertebral fracture, we test whether Tsimane fracture prevalence is higher than matched Americans, particularly for women. To determine whether Tsimane thoracic vertebral fracture results from compromised bone strength as opposed to trauma, we test whether fracture risk is inversely associated with thoracic vertebral BMD. Lastly, we test whether Tsimane women's fracture risk increases with reproductive effort, even after adjusting for BMD, which is expected if greater reproductive effort compromises bone micro-architecture in complex ways beyond just reducing mineral density (e.g. by reducing trabecular thickness or connectivity density).

Results
Tsimane women's thoracic vertebral BMD declines with early age at first birth and short IBI Earlier age at first birth is associated with reduced BMD (Std. b Age at 1st birth [years, logged] =0.099, p=0.036, controlling for age and fat-free mass, adj. R 2 = 0.51, n = 213; Appendix 1-table 1). Backtransforming logged age at first birth values into observed values and holding controls at sample means, there is a BMD difference of 0.57 SDs for women with maximum versus minimum age at first birth (37 versus 12 years, respectively).
Parity (continuously or categorically operationalized) is not associated with BMD controlling for age at first birth, age and fat-free mass (Appendix 1-table 2), nor does parity interact with any indicator of reproductive effort to predict BMD.
Shorter mean IBI (<29.7 months) is associated with lower BMD (Std. b Shorter mean IBI =-0.201, p=0.032, controlling for age at first birth, age and fat-free mass, adj. R 2 = 0.52; Appendix 1-tables 3-4). Mean IBI also interacts with age at first birth: BMD is 0.28 SDs higher for women with longer mean IBI and later age at first birth versus women with shorter mean IBI and earlier age at first birth (interaction p=0.027, controlling for age and fat-free mass; see Figure 1 and Appendix 1-figure 1).
Additionally controlling for indicators of modernization, that is, residential proximity to the closest market town of San Borja, Spanish fluency and schooling, which could reflect differential activity levels, diet and/or other factors affecting bone strength (e.g. infectious burden), strengthens the association between BMD and both mean IBI and age at first birth (comparing estimates in Appendix 1table 5 to those in Appendix 1-table 4). BMD is not significantly associated with any modernization indicator. Neither young age at menarche nor menopause is associated with BMD (Appendix 1-table 6), nor does either variable interact with any indicator of reproductive effort to predict BMD.
Lower thoracic vertebral BMD increases risk of thoracic vertebral fracture among Tsimane, particularly for women. Short IBI additionally increases Tsimane women's fracture risk Tsimane women with thoracic vertebral fracture (grade !1) have lower BMD (0.46 SDs; p=0.005), higher parity (0.25 SDs; p=0.071) and shorter mean IBI (0.22 SDs; p=0.011) than women without fracture (Appendix 1-table 12). There are no significant differences between women without versus with fracture in terms of age, anthropometrics, ages at menarche, menopause or first birth, or modernization indicators. Women's BMD is inversely associated with fracture risk (adjusted RR BMD per SD increase = 0.542, 95% CI: 0.352-0.837, p=0.006, controlling for age, height and fat mass, n = 219); this association remains (adjusted RR BMD = 0.540, 95% CI: 0.345-0.845, p=0.007; Figure 4) Figure 3. Thoracic vertebral (T6-T12) fracture prevalence (95% CIs) for Tsimane (in gray) and American (in green) women (A) and men (B) by fracture grade. Prevalence based on a less (grade !1) and more (grade !2) conservative fracture definition is shown. Prevalence is estimated from log-binomial generalized linear models adjusting for age. DOI: https://doi.org/10.7554/eLife.48607.004 after adding to the model mean IBI, which is also inversely associated with fracture risk (adjusted RR Mean IBI per SD increase = 0.379, 95% CI: 0.165-0.866, p=0.021; Appendix 1-table 13: Model 2). The inverse association between women's BMD and fracture risk strengthens with higher fracture grades (Appendix 1-table 14), whereas the inverse association between mean IBI and fracture risk is strongest for mild (i.e. grade 1) fracture. In multivariate models, we found no effect on women's fracture risk of either parity, ages at menarche, menopause or first birth, or modernization indicators (controlling for BMD, mean IBI, age, height and fat mass).
Tsimane men with fracture are shorter (0.23 SDs; p=0.038), have higher adiposity (0.37 SDs; p=0.025) and have higher BMI (0.34 SDs; p=0.012) than men without fracture (Appendix 1-table 15). As with women, men without versus with fracture do not significantly differ in terms of age, weight, fat-free mass or modernization indicators. But unlike women, there is no significant difference in BMD between men without versus with fracture, and men's BMD is not associated with risk of any fracture (adjusted RR BMD per SD increase = 1.046, 95% CI: 0.862-1.269, p=0.650, controlling for age, height and fat mass, n = 227; Appendix 1-table 16). Men's BMD is, however, inversely associated with risk of borderline (i.e. grade 0.5) deformity (Appendix 1-table 17; Figure 4). Modernization indicators are not associated with men's fracture risk, nor do they meaningfully moderate associations between BMD and fracture risk.

Discussion
This paper makes five empirical contributions. We find that: 1) Tsimane thoracic vertebral BMD is lower among women with early age at first birth and short IBIs; 2) Tsimane BMD is lower than a directly comparable American sample, but only for women; 3) thoracic vertebral fracture prevalence is higher for Tsimane than Americans; 4) among Tsimane, lower BMD is associated with higher fracture risk, particularly for women; and 5) short IBIs are associated with Tsimane women's fracture risk even after adjusting for BMD.
These results are consistent with a life history trade-off between reproductive effort and somatic maintenance (Kirkwood and Austad, 2000;Stearns, 1992;Stieglitz et al., 2015). Maternal physiology uniquely responds to the greater mineral demands of pregnancy and lactation by mobilizing . Association between thoracic vertebral BMD and thoracic vertebral (T6-T12) fracture risk (95% CI) for Tsimane. Log-binomial generalized linear models are used to estimate relative risk after adjustment for age, height and fat mass. Mean IBI is also included as a covariate for women. Parameter estimates are shown in Appendix 1-tables 13-14 (women), and Appendix 1-tables 16-17 (men). Severe fractures (grade 3) are omitted due to their relative scarcity. DOI: https://doi.org/10.7554/eLife.48607.005 skeletal mineral stores (Kovacs, 2016;Nelson et al., 2002;Prentice, 2003). Maternal regulatory mechanisms can compensate for acute mineral losses, for example, by retaining excess mineral in circulation to facilitate storage. But in energy-limited settings (due to high energy expenditure relative to consumption), earlier age at first birth and shorter birth spacing can potentially reduce peak bone mass and generate unbalanced cycles of bone resorption and formation, compromising bone strength long before menopause. In addition to direct negative effects of reproductive effort on bone strength, indirect negative effects are possible via reductions in women's PALs and mechanical loading on bone, given nursing women's reduced food acquisition efforts relative to non-nursing women (cf. Hurtado et al., 1992). Tsimane women's participation in certain subsistence activities ostensibly entailing high-impact and high-magnitude loadings (e.g. rice-pounding using thick wooden bats) may be curtailed during pregnancy and early lactation; meta analyses of controlled exercise trials indicate that high-impact and high-magnitude loadings are especially effective at preserving women's BMD at the lumbar spine, femoral neck and total hip (Martyn-St James and Carroll, 2009;Martyn-St James and Carroll, 2010). In the present study, we also found that parity per se is not inversely associated with BMD (Appendix 1-tables 2-3; but see Stieglitz et al., 2015). Prior studies of parity-specific effects on maternal BMD do not reveal consistent associations (Kovacs, 2016), although most studies are conducted in developed countries with lower fertility and greater energetic surpluses. These conditions can obscure expected energetic trade-offs between reproduction and maintenance.
The finding that thoracic vertebral BMD is lower for Tsimane than in a comparable sample of American women, whereas for men few population-level differences are apparent, despite higher mean PALs of Tsimane women and men, contributes to a growing literature emphasizing sexual dimorphism in skeletal responses to environmental stimuli (e.g. Macintosh et al., 2017), mediated in part by sex and growth hormones, growth factors and their receptors (Devlin, 2011;Gosman et al., 2011). Dietary or other systemic population-level differences (e.g. in infectious burden) do not easily explain Tsimane women's lower BMD because such systemic factors should affect both sexes. Recent dietary analyses of Tsimane and Americans (NHANES) based on 24 hr dietary recalls indicate minor sex differences in proportional macronutrient intake for both populations (Kraft et al., 2018). In absolute terms, Tsimane mean daily per-capita energy, protein and carbohydrate intake actually exceeds that of Americans for both sexes, as does Tsimane intake of the bone-forming minerals magnesium, phosphorus and zinc. Nevertheless, prior Tsimane research indicates that energetic limitation and greater immune activation from high pathogen exposure partly explain why mean PALs above those of industrialized societies do not inevitably yield elevated peak bone mass or protect against age-related bone loss for either sex (Stieglitz et al., 2015;Stieglitz et al., 2016). In a prior cross-sectional study of women utilizing ultrasound, we found reduced calcaneal strength for Tsimane versus Americans that is already apparent in the 20s, with population-level differences increasing with age (Stieglitz et al., 2015). Between-and within-population analyses in the present study suggest that sex-specific costs of reproduction contribute to skeletal sexual dimorphism, and that sexually dimorphic responses to reproduction manifest early in adulthood, depending on age at first birth and birth spacing (Madimenos et al., 2012;Stieglitz et al., 2015).
That Tsimane fracture prevalence is higher than a comparative American sample ( Figure 3 and Appendix 1-tables 8-11) may be surprising in light of the presumably protective higher lifetime moderate physical activity of Tsimane (Gurven et al., 2013), and their minimal exposure to other osteoporosis risk factors found in industrialized societies (e.g. glucocorticoid therapy, excessive smoking or alcohol consumption). By restricting the Tsimane sample to older adults (aged 40+ years), we minimize potential for cohort effects related to changing lifestyles associated with modernization. This is supported by the fact that study findings hold even after controlling for modernization indicators, which are not associated with thoracic vertebral BMD or fracture risk (Appendix 1-tables 5-6, 12, 15). Higher fracture prevalence among Tsimane versus Americans is noteworthy in light of relatively scant bio-archaeological evidence of osteoporotic fracture prior to industrialization (Agarwal, 2008;Curate et al., 2010) and evidence of increasing age-specific osteoporotic fracture incidence rates over time in Western populations (Cooper et al., 2011). But these results are not surprising in light of high Tsimane fertility and life history trade-offs between energetic investments in reproduction and bone growth and/or maintenance, and other factors (e.g. calcium deficiency, chronic immune activation due to high pathogen exposure), which may interact with high reproductive effort to further constrain the ability of bone tissue to respond to mechanical loading and high PALs (cf. Armelagos et al., 1972).
The fact that Tsimane women's vertebral BMD is inversely associated with vertebral fracture risk (cf. Mays, 1996) (Figure 4) suggests a major role of compromised bone strength in precipitating fracture, rather than traumatic injury of otherwise healthy bone. For Tsimane men, it is likely that trauma plays a major role in precipitating fracture due to: 1) a weak association between BMD and fracture risk ( Figure 4); 2) a high fracture prevalence relative to a comparative American sample (Figure 3) despite minimal BMD differences ( Figure 2); and our anecdotal observations of high levels of mechanical stress on Tsimane men's vertebrae from frequent heavy load carrying (e.g. of hunted game, timber for constructing houses). This of course does not preclude a contributing role of compromised bone strength (e.g. from micro-architectural deficiencies) in precipitating Tsimane men's fracture. Likewise, for women, the present results do not preclude a contributing role of trauma in precipitating fracture. Future research is needed to determine the extent to which women's subsistence involvement and frequent carrying of young children and other loads (e.g. woven bags filled with harvested cultigens) influence bone structural integrity. Some of the greatest compressive vertebral loads occur when weights are carried in front of the body (Bruno et al., 2017;Rohlmannt et al., 2001), which is how Tsimane women routinely carry infants and toddlers. The fact that Tsimane women's shorter IBIs predict increased fracture risk even after adjusting for BMD, which also remains a significant predictor (Appendix 1-table 13-14), suggests that shorter IBIs compromise multiple aspects of vertebral micro-architecture (e.g. trabecular thickness or connectivity density), although further research is needed to examine this possibility.

Inferring behaviors underlying morphological variation in past human populations
While our results do not directly address debates over the timing of and prior selection pressures underlying transition to skeletal gracility in past human populations, our results do provide insight into lifestyle factors affecting bone strength which may have been relevant during this transition. It has been hypothesized that subsistence transition from hunting and gathering to more sedentary agriculture, and increasing reliance on labor-saving technology, caused reductions in mechanical loading on bone and PALs, leading to modern human skeletal gracilization (Chirchir, 2019;Chirchir et al., 2015;Chirchir et al., 2017;Ruff et al., 2015;Ryan and Shaw, 2015). Increasing reliance on agriculture may have been associated with reductions in terrestrial mobility, and thus reduced mechanical loading of lower limbs, although the pace and magnitude of these changes likely varied temporally and spatially. Changes in upper body activities may have been much more variable during subsistence transitions, so upper limb loading may have actually increased with agriculture in some regions (Bridges, 1989). Nevertheless, comparisons of skeletal remains of huntergatherers and either full-or part-time agriculturalists indicate among agriculturalists reduced femoral strength, as indicated by trabecular bone structure or external size dimensions (Larsen, 1981;Ryan and Shaw, 2015), and accelerated age-related decline in radial bone mineral content (Perzigian, 1973).
Yet a physical inactivity explanation for modern human skeletal gracility -rooted in a subsistence transition from foraging to farming -is puzzling for several reasons. Evidence that agriculturalists are more sedentary than hunter-gatherers is not particularly strong: PALs and time allocation to work vary substantially within a subsistence regime, and both measures are actually higher among agriculturalists (Gurven et al., 2013;Leonard, 2008). Moreover, children in agricultural societies generally begin work earlier than hunter-gatherer children (Kramer, 2005). This is significant because higher PALs in childhood and early adulthood, particularly for higher-impact activities producing high peak stresses, increase peak bone mass, size and later-life bone strength (e.g. see Warden et al., 2014 and references therein). Furthermore, in a prospective Tsimane study we found that time spent in horticulture positively predicted ultrasound-derived indicators of radial strength , whereas tibial strength was not predicted by extent of involvement in any subsistence activity including hunting. While various studies show reduced lower limb strength among agriculturalists relative to pre-agriculturalists, as indicated by cross-sectional diaphyseal structure of cortical bone or trabecular bone volume fraction or thickness (May and Ruff, 2016;Saers et al., 2016), other studies using similar diaphyseal structural properties or external size dimensions show no differences in limb strength or dissimilar patterns by subsistence regime (Bridges, 1989;Ruff, 1999). Taken together, these observations create uncertainty over the timing of and selection pressures underlying modern human skeletal gracility. This uncertainty is exacerbated given evidence of sex-specific skeletal responses to mechanical loading and physical activity (Bridges, 1989;Macintosh et al., 2017;Ruff, 1999). Additional sex-specific factors may have contributed to skeletal gracility during transition to farming, either independently and/or in interaction with PALs.
Our results instead suggest that fertility increases associated with subsistence transition from foraging to farming (Bentley et al., 1993) contributed to modern human skeletal gracility, particularly for women. Numerous proximate determinants have been proposed to explain fertility increases that accompanied greater energetic surpluses from agriculture, including earlier age of menarche and first birth, and shorter IBIs (e.g. due to earlier weaning and supplementary infant feeding) (e.g. Campbell and Wood, 1988). Regardless of proximate fertility determinants, given the trade-off between energetic investment in reproduction and somatic maintenance, our results suggest that women's skeletons were especially susceptible to gracilization during subsistence transition from foraging to farming in light of the associated increases in fertility, and reduced mobility resulting from increased fertility (cf. Hurtado et al., 1992). Sex-specific mechanisms beyond menopause underlying skeletal gracility during this subsistence transition should thus be considered, in addition to explanations emphasizing increased sedentism, reliance on labor-saving technology and associated reductions in mechanical loading on bone. Nevertheless, since skeletal gracilization is observed in both sexes during this subsistence transition , any explanation solely resting on changes in reproductive effort cannot explain these morphological changes in men.

Study limitations
The cross-sectional study design using retrospective demographic data limits our ability to establish that greater reproductive effort causes BMD reductions and fracture. The results presented here may thus be consistent with alternative interpretations derived from life history theory or proximate explanations, although such alternatives must address the observed population-level sex differences. Another limitation is that our measures of reproductive effort are indirect measures of reproductive costs, and we lack data on lactation duration or intensity. However, all Tsimane women breastfeed their infants and it is common for women to breastfeed exclusively for about four months (no study participant bottle fed an infant). We also lack estimates of vertebral size and geometry, which affect vertebral strength, and we lack individual-level data on activity level, nutrient intake and pathogen burden for population-level comparisons of BMD and fracture prevalence and severity. We also do not consider whether genetic diversity determines heterogeneity in BMD or fracture prevalence. Bone strength indicators such as BMD and fracture are heritable and the frequency of alleles affecting BMD differs between ethnically distinct populations (see Wallace et al., 2016 and references therein). But while genetic factors can account for a sizable portion of variance in bone strength within populations, there is little evidence that heterogeneity in bone strength between populations is due to stochastic genetic diversity.

Conclusion
This study examines direct indicators of bone strength from a clinically and mechanically relevant anatomic region using in vivo imaging in a physically active population with high fertility. Results suggest a trade-off between reproductive effort and bone strength, and that greater reproductive effort constrains the ability of bone tissue to respond to mechanical loading and high physical activity. Results also raise the possibility that increased fertility associated with subsistence transition from foraging to farming promoted modern human skeletal gracility, particularly among women. Because of the complex nature of lifestyle transformations during subsistence transitions, including apparent increases in infectious disease, nutritional deficiencies and dental decay (Cohen and Armelagos, 1984), an expanded conceptual framework incorporating diverse lifestyle factors that may constrain the ability of bone to respond to mechanical loading (e.g. high fertility, nutrient deficiency, infectioninduced inflammation (Madimenos et al., 2012;Stieglitz et al., 2015;Stieglitz et al., 2016) can improve our understanding of morphological transformations associated with transition to farming. Of course, our ability to make inferences about the past using data collected in contemporary populations is limited. Tsimane are neither 'pure' hunter-gatherers nor agriculturalists, and they may differ in important ways from ancestral human populations in terms of residential mobility, fertility, diet and disease exposures. Yet no single population represents the range of experiences across different environments that shaped the evolution of our species over the millennia in which ecologies fluctuated. In vivo study of bone strength in well-characterized, population-representative, nonindustrialized societies provides an opportunity to examine lifestyle factors that are often invisible to bio-archaeological inquiry but nonetheless relevant to understanding selection pressures over human history.

Study design and participants
The Tsimane sample includes all individuals who met the inclusion criteria of self-identifying as Tsimane and who were aged 40+ years (n = 507; 48% female; age range: 41-94 years; see Appendix for additional details and Appendix 1-table 18 for descriptives of all study variables). 185 of the 245 participating women (76%) were post-menopausal. No participant reported ever using hormonal contraception or dietary supplements with consistency. No Tsimane was excluded based on any health condition that can affect BMD or fracture risk.
Comparative American BMD data were collected among asymptomatic subjects from greater Los Angeles as part of a different study (described in Budoff et al., 2010). Briefly, 9585 subjects (43% female; mean age = 56) underwent coronary artery calcification (CAC) scanning for evaluation of subclinical atherosclerosis, after exclusion of participants with vertebral deformities or fractures. Subjects had no known bone disease (see Appendix 1-table 7 for additional details). Two American data sources are used to compare Tsimane and American fracture prevalence: a subset from the MESA study (Budoff et al., 2011) and a subset reported in Budoff et al. (2013) (see Appendix 1table 8 for additional details); data from these American subsets were matched to Tsimane by age, sex and weight (±5 kg), and then merged to create a single American comparison sample.

Thoracic computed tomography (CT)
Tsimane CT scans were conducted at the Hospital Presidente German Busch in Trinidad, Bolivia using a 16-detector row scanner (GE Brightspeed, Milwaukee, WI, USA). A licensed radiology technician acquired a single, ECG-gated non-contrast thoracic scan as part of a broader project on atherosclerosis, including CAC assessment (see  and Appendix for additional details). Typical multi-detector CT protocols used for evaluation of CAC include imaging the mid-thoracic spine in the reconstructed field of view, facilitating simultaneous evaluation of thoracic vertebrae during a single examination without additional radiation exposure. Tsimane CT settings were: 250 ms exposure, 2.5 mm slice thickness, 0.5 s rotation speed, 120 kVp, and 40 mA with prospective triggering. Refer to the Appendix for details on CT parameters for the comparative American samples; American CT data were all collected at the same institute (Los Angeles Biomedical Research Institute).
Vertebral BMD was measured manually in each of three consecutive thoracic vertebrae (T7-T10 range) by a radiologist with 20+ years of experience (see Appendix for additional details). BMD measurement started at the level of the section that contained the left main coronary artery (LMCA) caudally (beginning at either T7 or T8, depending on the origin of the LMCA). The center of the region of interest was located at the center of each vertebrae, with a 2-3 mm distance from the cortical shell; this distance ensured that BMD measurements within the vertebral body excluded the cortical bone of the vertebral shell. For each vertebrae, the radiologist manually positioned a circular region of interest while demarcating cortical from trabecular bone based on visual inspection. Any area with large vessels, bone island fractures and calcified herniated disks were excluded as much as possible from the region of interest with use of the manual free tracing protocol. Mean BMD for the three consecutive thoracic vertebrae was then calculated. This BMD measure is strongly positively correlated (Pearson r's > 0.9) with lumbar vertebral BMD (Budoff et al., 2012). CT-derived BMD estimates can be obtained with and without calibration phantoms. Phantomless BMD estimates correlate strongly (Pearson r = 0.99) with standard phantom-based CT BMD estimates (Budoff et al., 2013). Hounsfield units were converted to BMD (mg/cm 3 ) using a calibration phantom of known density or a scanner-specific mean calibration factor for the T7-T10 vertebrae from scans performed without the phantom. All BMD measurements used in this study were performed at the Los Angeles Biomedical Research Institute.

Thoracic vertebral fracture
For each subject the radiologist classified seven vertebrae (T6-T12) according to Genant's semiquantitative technique (GST) (Genant et al., 1993). While there is no consensus regarding the radiologic definition of vertebral fracture, the GST provides highly reproducible diagnosis of fractures, is the current clinical technique of choice for diagnosing fracture, and is the most widely used technique for identifying fracture (Shepherd et al., 2015). Based on visual inspection, each vertebra is rated according to severity of loss of vertebral height and other qualitative features, including alterations in shape and configuration of the vertebra relative to adjacent vertebrae and expected normal appearances. Each vertebra is classified into one of five categories: normal (grade 0); mild fracture (grade 1; approximately 20-25% reduction in anterior, middle, and/or posterior vertebral height, and a 10-20% reduction in projected vertebral area); moderate fracture (grade 2; 25-40% reduction in any height and a 20-40% reduction in area); and severe fracture (grade 3; >40% reduction in any height and area). A grade 0.5 indicates borderline deformed vertebra (<20% reduction in any height) that is not considered to be a definitive fracture (see Appendix for additional details). Each subject is assigned one grade representing a summary measure of all seven vertebrae. Subjects with >1 vertebral deformity are classified according to their most severe deformity. Subjects are considered to present vertebral fracture if any vertebral body is graded at least mildly deformed (i.e. grade !1); subjects are considered to present no fracture if graded 0 or 0.5. Given recent analyses (Lentle et al., 2018) showing lower observer agreement for mild fractures (the most common) relative to moderate and severe fractures, we repeat analyses using a more conservative fracture definition (i.e. grade !2). All fracture measurements used in this study were performed at the Los Angeles Biomedical Research Institute.

Socio-demographics and anthropometrics
Individuals for whom reliable ages could not be ascertained are not included in analyses (see Appendix for additional details). Reproductive histories were elicited in the Tsimane language. IBI refers to the number of months between live births for women with !2 live births. Self-reported ages at menarche and menopause were recorded during medical exams conducted by physicians of the Tsimane Health and Life History Project (THLHP). During annual THLHP census updates we also coded for each participant their village of residence (from which we derived via GPS residential proximity to the closest market town of San Borja), self-reported Spanish fluency (0 = none; 1 = moderate; 2 = fluent) and schooling (# years) as indicators of modernization.
Height and weight were measured during THLHP medical exams using a Seca stadiometer (Road Rod) and Tanita scale (BC-1500). The scale uses a method of bioelectrical impedance analysis to estimate percent body fat. Using weight and percent body fat we calculated fat mass (weight*percent body fat) and fat-free mass (weight -fat mass).

Data analysis
The two outcome variables indicating bone strength are thoracic vertebral BMD and thoracic vertebral fracture. Sexes are analyzed separately given the sex-specific nature of hypotheses and to minimize confounding by unobserved factors. General linear models are used to test for associations between Tsimane women's BMD and reproductive effort (see Appendix for additional details). We compared BMD of Tsimane and age-and sex-matched Americans using age-standardized means. We use a parametric test (one sample t test) to evaluate whether population-level differences in mean BMD within each decade are significant (p<0.05), specifying as the test value the US means published in Budoff et al. (2010). We compare Tsimane and American fracture prevalence using age-standardized values, and using log-binomial generalized linear models (GLMs). GLMs are used to test for effects of BMD and women's reproductive effort on the probability of fracture. Both continuous and categorical (e.g. median split) measures of reproductive effort are used to analyze their associations with BMD and fracture risk. Unless otherwise noted, anthropometric and socio-demographic covariates are included in a stepwise fashion in regressions (see Appendix for descriptive analyses of Tsimane BMD by age and sex, and by anthropometrics [Appendix 1-tables 19-20]). Participants with any missing values are removed from analyses. The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.

Data availability
The data that support the findings of this study are available on Dryad (https://doi.org/10.5061/ dryad.rf0g0md).
The following dataset was generated: Author (    a US data are from two sources: 1) a subset of MESA study participants (described in Budoff et al., 2011), and 2) a subset of study participants in the greater Los Angeles area (described in Budoff et al., 2013). were 1214 eligible people living in these villages (the only eligibility criteria were being 40+ years old, self-identifying as Tsimane and willing to participate). 731 individuals were present in their villages at the time and subsequently received a CT scan. Transporting participants from their village to the nearby market town of San Borja was logistically complicated (requiring trekking through the forest, dug-out canoes, rafts propelled by poles pushed off the river bottom, trucks, and cars) and can require up to two days of travel each way. From San Borja to the Beni department capital of Trinidad (where the hospital containing the CT scanner is located) is an additional 6 hr car ride. Due to these logistical complications, participants not in their village at the time we arrived were not sampled. The Tsimane are semi-mobile and often build secondary houses deep in the forest near their horticultural plots, not returning to their village for extended periods of time. Hunting and fishing trips can last days or weeks, and some men engage in wage labor in San Borja or elsewhere (e.g. rural cattle ranches). In an average village, approximately one-third of individuals are away hunting, fishing, working in their horticultural plots, or in San Borja at any given time.
Additionally, a major flood in February 2014 resulted in mass migration from some villages, and the creation of several new villages that were not sampled as part of this study, further reducing the number individuals that could be sampled in this study. To address potential sources of sample bias, analyses comparing Tsimane who received CTs and those who did not but who participated in the THLHP's medical exams in Tsimane villages were conducted. There were no significant differences in sex, blood pressure or body fat (see Kaplan et al., 2017 for additional participant details) and thus the CT sample analyzed here is thought to be representative of all Tsimane aged 40+ years.
Of the 731 Tsimane who received a CT scan, CT data from 507 (69%) were selected with no particular criteria to estimate thoracic vertebral bone mineral density (BMD) and fracture (scans from 224 Tsimane were not assessed due to radiologist time constraints). Among these 507 individuals, some data were missing because of either broken equipment, missing supplies, participant recall problems, absent or sick team personnel who were unable to collect data, or because all relevant thoracic vertebrae did not appear in the CT image field of view (see Appendix 1- Thoracic computed tomography (CT) Tsimane CT scans were supervised and reviewed by at least one of the HORUS study team cardiologists and radiologists. Breath-hold instructions were given in the Tsimane language to minimize respiratory motion artifact and misregistration (for additional methodological details see . CT scanning procedures used to collect BMD data from greater Los Angeles involved use of electron-beam CT scanners (C-300; GE-Imatron, South San Francisco, California) and a 64-detector row CT scanner (LightSpeed VCT; GE Medical Systems) at the Los Angeles Biomedical Research Institute. Parameters for electron-beam CT scanning were 130 kVp, 630 mA, and 2.5 or 3 mm slice thickness. The multidetector CT parameters were 120 kVp, 200-600 mA, and 2.5 mm slice thickness. 156 subjects underwent scanning with electron-beam CT and 64-detector row CT on the same day for comparison and normalization of BMDs between scanners. CT scanning to collect MESA sub-sample data, used here to assess fracture prevalence and severity, was performed at the Los Angeles Biomedical Research Institute; procedures involved the use of a C-150 GE-Imatron electron-beam scanner and were otherwise identical to those described in Budoff et al. (2010). Supplementary CT data for assessment of fracture prevalence and severity were also collected at the Los Angeles Biomedical Research Institute using similar parameters.

Thoracic vertebral bone mineral density (BMD)
CT-measured vertebral BMD is increasingly used for osteoporosis screening because of its ability to provide three-dimensional information compared to traditional dual x-ray absorptiometry two-dimensional images (Budoff et al., 2012). The left main coronary artery (LMCA) was set as the reference site to allow reproducible detection of a spinal level for use with cardiac CT scanning. The LMCA is covered in 100% of images and the field of view can be completely reconstructed; therefore, it is an optimal reference point to locate the starting measurement level. The most common location of the LMCA origin is at the lower edge of T7.
Intra-observer variability in BMD measurements using the same measurement technique in a different sample (Budoff et al., 2010) was tested on 120 scans by one observer from the Los Angeles Biomedical Research Institute, with 1 week intervals between the two readings. To measure inter-observer variability on 67 randomly selected scans, the results obtained by two radiologists from the Los Angeles Biomedical Research Institute who were blinded to all clinical information and prior measurements were compared using this other sample. Intraand inter-observer variations in BMD measurements were 2.5% (Bland-Altman plot ratio: 1.00, 95% CI: 0.99-1.00) and 2.6% (Bland-Altman plot ratio: 1.00, 95% CI: 0.99-1.00), respectively.

Thoracic vertebral fracture
Genant's semi-quantitative technique does not distinguish between wedge (i.e. reduced anterior height), biconcave (i.e. reduced central height) or crush (i.e. reduced posterior height) fractures; most fractures contain a combination of these features and are influenced by the local biomechanics of the spinal level involved. Vertebral fractures are differentiated from other, non-fracture deformities (e.g. osteoarthritis), although these other deformities were not systematically coded.

Socio-demographics and anthropometrics
Birth years were assigned based on a combination of methods described in detail elsewhere (Gurven et al., 2007), including using known ages from written records, relative age lists, dated events, photo comparisons of people with known ages, and cross-validation of information from independent interviews of kin. Each method provides an independent estimate of age, and when estimates yielded a date of birth within a three-year range, the average was generally used.
Outcomes of each pregnancy reported during reproductive histories were recorded as either ending in a live birth or terminating pre-term. Whether miscarriages (including stillbirths) are included or omitted from parity counts does not affect results, and results reported here reflect only live births.

Data analysis
With respect to general linear models of BMD, variance inflation factors were checked to assess degree of covariance and ensure minimal collinearity. Crude mean BMD values are used to compare populations for each decade, and age-standardized means are used to compare populations across all decades. We use the Tsimane adult age distribution, estimated from the THLHP 2015 census, as the standard.
Appendix 1-table 19. General linear models: anthropometric predictors of thoracic vertebral BMD (mg/cm 3 ) for women. Height is included in models 1 and 3, weight in models 2, 3 and 5, body mass index (BMI) in model 4, body fat (% or fat mass) in models 5-6, and fat-free mass in model 6. Age is controlled in all models. Standardized betas are shown (intercepts omitted).