An integrative skeletal and paleogenomic analysis of stature variation suggests relatively reduced health for early European farmers.

Human culture, biology, and health were shaped dramatically by the onset of agriculture ∼12,000 y B.P. This shift is hypothesized to have resulted in increased individual fitness and population growth as evidenced by archaeological and population genomic data alongside a decline in physiological health as inferred from skeletal remains. Here, we consider osteological and ancient DNA data from the same prehistoric individuals to study human stature variation as a proxy for health across a transition to agriculture. Specifically, we compared “predicted” genetic contributions to height from paleogenomic data and “achieved” adult osteological height estimated from long bone measurements for 167 individuals across Europe spanning the Upper Paleolithic to Iron Age (∼38,000 to 2,400 B.P.). We found that individuals from the Neolithic were shorter than expected (given their individual polygenic height scores) by an average of −3.82 cm relative to individuals from the Upper Paleolithic and Mesolithic (P = 0.040) and −2.21 cm shorter relative to post-Neolithic individuals (P = 0.068), with osteological vs. expected stature steadily increasing across the Copper (+1.95 cm relative to the Neolithic), Bronze (+2.70 cm), and Iron (+3.27 cm) Ages. These results were attenuated when we additionally accounted for genome-wide genetic ancestry variation: for example, with Neolithic individuals −2.82 cm shorter than expected on average relative to pre-Neolithic individuals (P = 0.120). We also incorporated observations of paleopathological indicators of nonspecific stress that can persist from childhood to adulthood in skeletal remains into our model. Overall, our work highlights the potential of integrating disparate datasets to explore proxies of health in prehistory.

The agricultural revolution—beginning ∼12,000 B.P. in the Fertile Crescent zone (1, 2) and then spreading (3–5) or occurring independently (6, 7) across much of the inhabited planet—precipitated profound changes to human subsistence, social systems, and health. Seemingly paradoxically, the agricultural transition may have presented conflicting biological benefits and costs for early farming communities (8, 9). Specifically, demographic reconstructions from archaeological and population genetic records suggest that the agricultural transition led to increased individual fitness and population growth (6, 10–12), likely due in part to new food production and storage capabilities. Yet, bioarchaeological analyses of human skeletal remains from this cultural period suggest simultaneous declines in individual physiological well-being and health, putatively from 1) nutritional deficiency and/or 2) increased pathogen loads as a function of greater human population densities, sedentary lifestyles, and proximity to livestock (9, 13–18).

To date, anthropologists have used two principal approaches to study health across the foraging-to-farming transition in diverse global regions (13, 19, 20). The first approach involves identifying paleopathological indicators of childhood stress that persist into adult skeletal remains. For example, porotic hyperostosis (porous lesions on the cranial vault) and cribra orbitalia (porosity on the orbital roof) reflect a history of bone marrow hypertrophy or hyperplasia resulting from one or more periods of infection, metabolic deficiencies, malnutrition, and/or chronic disease (21–26). Meanwhile, linear enamel hypoplasia (transverse areas of reduced enamel thickness on teeth) occurs in response to similar childhood physiological stressors (e.g., disease, metabolic deficiencies, malnutrition, weaning) that disrupt enamel formation in the developing permanent dentition (27–30). Broadly, these paleopathological indicators of childhood stress tend to be observed at higher rates among individuals from initial farming communities relative to earlier periods, potentially reflecting their overall “poorer” health (14, 31–36).

A second approach uses skeleton-based estimates of achieved adult stature as a proxy for health during childhood growth and development (37–39). Since stature is responsive to the influences of nutrition and disease burden alongside other factors, relatively short “height-for-age” (or “stunting”) has been used as an indicator of poorer health in both living and bioarchaeological contexts (39–43). When studying the past, individual stature can be estimated from long bone measurements and regression equations (44–47). Using these methods, multiple prior studies have reported a general profile of relatively reduced stature for individuals from early agricultural societies in Europe (15, 48–50), North America (51–53), the Levant (16, 32), and Asia (54, 55). For example, estimated average adult mean statures for early farmers are ∼10 cm shorter relative to those for preceding hunter-gatherers in both western Europe (females, −8 cm; males, −14 cm) (49, 50) and the eastern Mediterranean (females, −11 cm; males, −8 cm) (56). This pattern is not universal, as a few studies do not report such changes (57, 58); the variation could be informative with respect to identifying potential underlying factors (59).

However, in addition to environmental effects like childhood nutrition and disease, inherited genetic variation can have an outsized impact on terminal stature, with ∼80% of the considerable degree of height variation within many modern populations explainable by heritable genetic variation (60–63). Moreover, migration and gene flow likely accompanied many subsistence shifts in human prehistory. For example, there is now substantial paleogenomic evidence of extensive population turnover across prehistoric Europe (64–69). Therefore, from osteological studies alone, we are unable to quantify the extent to which temporal changes in height reflect variation in childhood health vs. changes/differences in the frequencies of alleles associated with height variation.

In this study, we have performed a combined analysis of ancient human paleogenomic and osteological data where both are available from the same n = 167 prehistoric European individuals representing cultural periods from the Upper Paleolithic (∼38,000 B.P.) to the Iron Age (∼2,400 B.P.). This approach allows us to explore whether “health,” as inferred from the per-individual difference between predicted genetic contributions to height and osteological estimates of achieved adult height, changed over the Neolithic cultural shift to agriculture in Europe. When craniodental elements were preserved and available for analysis (n = 98 of the 167 individuals), we also collected porotic hyperostosis, cribra orbitalia, and linear enamel hypoplasia paleopathological data in order to examine whether patterns of variation between osteological height and genetic contributions to height are explained in part by the presence/absence of these indicators of childhood or childhood-inclusive stress.

Results

We developed a database of n = 167 ancient European adult individuals (67 females, 100 males) with available genome-wide paleogenomic data (from either shotgun sequencing or DNA capture-based approaches and from both published and in-process studies) and stature estimates based on long bone measurements (either newly collected or published) (Fig. 1 and Dataset S1). The cultural and time periods represented in our dataset are Upper Paleolithic (38,000 to 12,000 B.P.; n = 8 individuals), Mesolithic (11,000 to 6,400 B.P.; n = 15), Neolithic (7,100 to 3,500 B.P.; n = 46), Copper Age (6,300 to 3,400 B.P.; n = 60), Bronze Age (4,500 to 2,500 B.P.; n = 31), and Iron Age (2,600 to 2,400 B.P.; n = 7).

Fig. 1.

Osteological stature and ancient DNA–based polygenic height scores. (A) Map of the locations of the archaeological sites from which individuals included in the dataset were recovered. (B) Osteological height estimates generated using measurements of long bone lengths (highlighted in red on the illustration) and sex-specific regression equations (44). (C) Polygenic height scores generated using genome-wide association summary statistics for height-associated single-nucleotide polymorphisms and individual ancient DNA genotype data. (D and E) The relationship between polygenic height score and estimated osteological stature (centimeters) for females, for males, and for the full sample with height differences from mean stature calculated separately for females (mean = 154.64 ± 6.48 cm) and males (mean = 165.97 ± 6.60 cm), respectively (represented on the y axis in E). (F) Residuals of the relationship between polygenic height score and osteological height with sex as a covariate for all individuals by cultural period. Mean and median are represented by the black and blue dashed lines, respectively. Skeletal illustration in A image credit: Katharine Thompson (Stony Brook University, Stony Brook, NY).

Subsistence strategies in the Upper Paleolithic and Mesolithic focused on gathering, collecting, and hunting food. The Neolithic is marked by the emergence of plant cultivation and animal domestication (to varying degrees and tempos of integration), long-term settlements, larger populations, and increased social complexity—processes that then intensified and expanded in subsequent periods (70). The overlapping dates among the different cultural periods reflect both geographical variation in the timing of cultural change and the potential for co-occurrence of multiple cultural traditions within a single region.

Confirmation of an Average Osteological Stature Dip in the Neolithic.

We used an osteometric board to newly estimate the lengths of preserved long bones for n = 93 of the n = 167 total individuals (55.7%) in our database (Dataset S1). We also recorded published and unpublished (previously collected) long bone length estimates for an additional n = 54 individuals (32.3%). In these cases, we estimated osteological stature from the long bone length data (44). Finally, for n = 20 individuals, only precalculated terminal height estimates were available (12%).

We observed osteological stature variation among cultural periods (Fig. 1). Reconstructions of osteological stature for females and males are lower during the initial shift to farming during the Neolithic compared with earlier and post-Neolithic periods. Specifically, individuals from the pre-Neolithic periods (female average stature = 157.28 ± 7.0 (SD) cm, males = 168.6 ± 7.6 cm) were ∼4 cm taller on average than those from the Neolithic (females = 151.0 ± 6.7 cm; males = 164.5 ± 7.4 cm; linear model including sex as a factor; P = 0.012). Then, Neolithic individuals were ∼2 to 5 cm shorter on average compared with all post-Neolithic individuals (females = 155.06 ± 5.7 cm; males = 166.16 ± 5.82 cm; P = 0.046), with average osteological stature steadily rebounding across the Copper Age (females = 154.8 ± 6.2 cm; males = 165.5 ± 5.5 cm), Bronze Age (females = 155.6 ± 4.8 cm; males = 167.3 ± 6.3 cm), and Iron Age (females = 158.2 ± 7.8 cm; males = 165.8 ± 1.3 cm) ().

The overall pattern from our data roughly parallels previously published reports. Specifically, 1) stature decreased slightly from the Upper Paleolithic to the Mesolithic (35, 36, 71), 2) marked stature reduction occurred during the initial agricultural transition in the Neolithic (14, 16, 34, 55) [although this is not universal (70, 72–74)], and 3) stature rebounded during subsequent post-Neolithic periods of agricultural intensification (14, 16).

Early Farmers Were Relatively Shorter than Expected Given Their Polygenic Height Scores.

We next considered the osteological height estimates in the context of ancient DNA-based polygenic height scores for the same individuals. Using an established approach for working with ancient DNA genotype data (65, 75, 76), we estimated a polygenic score for each prehistoric individual based on their available genome-wide genotypes in the context of results from a large-scale genome-wide association study (GWAS) of stature variation in modern Europeans (77) (data are from http://www.nealelab.is/uk-biobank/). For the results presented in the text and figures, we used a version of the dataset in which all variants possibly affected by deamination-based ancient DNA damage (78–80) were masked. We also performed all analyses with the unmasked dataset and obtained similar results ().

While the polygenic scores that we estimated for the n = 167 individuals were somewhat variable across cultural periods (Fig. 1 and ), as anticipated based on prior work (58), we were most interested in using these data to begin to account for genetic contributions to achieved adult (osteological) height on a per-individual basis. Polygenic height scores and osteological estimates of stature were positively correlated for females (n = 67 [Fig. 1; r2 = 0.076; P = 0.014), for males (n = 100 [Fig. 1; r2 = 0.073; P = 0.004), and for the combined dataset (Fig. 1) (r2 = 0.0796; P = 0.0001). These results support the general biological plausibility of our integrative analysis of paleogenomic and osteological data.

Importantly and expectedly, there is still considerable interindividual variation in the relationship between polygenic height score and achieved adult stature, which could reflect any combination of incomplete genetic information; long bone measurement error; polygenic height score or stature estimate error; and the effects of childhood nutrition, disease, and other environmental variables on growth. Accordingly, we next analyzed the residuals from the combined sex osteological stature and ancient DNA-based polygenic score model (Fig. 1) to test whether individuals tended to have taller or shorter adult stature relative to expectations given their individual polygenic scores across the different cultural periods. These residuals are expressed in plus or minus centimeters from the predicted stature per individual.

We present our results as differences between average residuals (reflecting the difference between osteological stature and genetic height score after accounting for sex) across cultural periods. As such, our below use of the term “expected” in relative statements in this context reflects a statistical property of the overall dataset rather than commentary on the broader anthropological hypotheses being tested.

When using our ancient DNA-based approach to partly account for the predicted contribution of genetic variation to adult stature, we observed that individuals from the Neolithic were indeed osteologically shorter than expected (i.e., based on their polygenic height scores and in the context of our overall sample) compared with individuals from other cultural periods (Fig. 1 and ). Specifically, pre-Neolithic individuals (average residual = +1.96 ± 7.06 cm) were +3.82 cm taller than expected on average relative to Neolithic individuals (average residual = −1.87 ± 7.08 cm; P = 0.040). Neolithic individuals were then −2.21 cm shorter than expected on average relative to post-Neolithic individuals (average residual = 0.341 ± 5.48 cm; P = 0.068). Neolithic individuals were the only group with a negative average residual.

We confirmed that these results cannot be explained by geographic variation (latitude and longitude) in our sample (). We also obtained similar results when we separately analyzed females and males () and when we separately analyzed the lengths of individual long bones as opposed to the reconstructed stature estimates (for the femur and radius, in particular, although there are sample size limitations) ().

In contrast, our results were partially muted when we included variables explicitly reflecting genetic ancestries in our model. Our primary approach (as above) already accounts for ancestry variation via individual-level calculations of polygenic scores. However, polygenic height scores explain only a proportion of total heritable variation (81–83). Therefore, we repeated our hypothesis tests after conducting a multidimensional scaling (MDS) analysis with the genome-wide genotype data for all n = 167 ancient individuals to then including the first four MDS components from this analysis as factors in an updated linear model, following Cox et al. (84). Similar to principal components analysis (PCA), MDS can be used to mitigate the effects of ancestry heterogeneity in a GWAS framework (85, 86). When including the MDS components in the model, pre-Neolithic individuals (average residual = +2.22 ± 6.98 cm) were +2.82 cm taller than expected on average relative to Neolithic individuals (average residual = −0.594 ± 6.9 cm; P = 0.12). Neolithic individuals were 0.38 cm shorter on average relative to post-Neolithic individuals (average residual = −0.21 ± 5.16 cm; P = 0.74) ().

Paleopathological Indicators of Nonspecific Stress.

Adverse early life conditions may negatively impact adult stature. To begin to investigate whether individual-level early life effects on prehistoric stature could be identified, we incorporated observations of paleopathological indicators of nonspecific stress that can persist from childhood to adult skeletal remains into our analytical model. To do so, we characterized the presence/absence of one or more of cribra orbitalia (porosity on the orbital roof), porotic hyperostosis (porosity on the cranial vault), and linear enamel hypoplasia (reduced areas of enamel thickness) for n = 98 of the n = 167 (58.7%) individuals in our study (n = 82 newly characterized; n = 16 published/previously characterized) ().

For 58 of these 98 individuals (59.2%), crania were sufficiently complete for assessment of the presence/absence of all three stress indicators (Fig. 2). Of this subsample, one or more stress indicators were present for 41 (70.7%) of the individuals, two or more indicators were observed in 18 (31.0%) individuals, and all three paleopathological indicators were present in only 2 (3.4%) individuals. Thus, stresses on health were relatively common overall in prehistoric Europe.

Fig. 2.

Paleopathological indicators of stress. Paleopathological indicators of nonspecific stress evaluated in this study: linear enamel hypoplasia (LEH; bands of reduced enamel thickness on teeth), cribra orbitalia (CO; porosity on the orbits), and porotic hyperostosis (PO; porosity on the side of the skull). (A) The remains of 58 individuals were sufficiently complete to permit presence/absence assessment for all three paleopathologies. The proportions of individuals with one or more, two or more, and all three stress indicators are indicated across cultural periods. Numbers above the bars indicate sample sizes. (B) The presence/absence of at least one of the three paleopathological indicators of stress could be determined for 98 total individuals. Shown are the proportions of individuals (of those who could be assessed for that indicator) with LEH, CO, and PH across cultural periods. Numbers above the bars indicate sample sizes. (C) Residuals of the relationship between polygenetic height score and osteological height with sex as a covariate plotted separately for individuals with each paleopathological indicator of stress present vs. absent. Means are represented by the black lines. Numbers above the bars indicate sample sizes. (D) Residuals of the relationship between polygenetic height score and osteological height with sex as a covariate plotted separately for individuals with and without cribra orbitalia by cultural period. Means are represented by the black lines. Numbers above the bars indicate sample sizes. Skeletal illustrations image credit: Katharine Thompson (Stony Brook University, Stony Brook, NY).

A considerable 77.8% (14 of 18) of Neolithic individuals had one or more stress indicators (Fig. 2). While the proportion of Copper Age individuals with one or more stress indicators (10 of 18; 55.6%) was lower compared with that for the Neolithic (Fisher’s exact test; P = 0.289), the Neolithic result is not unique, with one or more stress indicators also recorded for all but one individual in the Bronze Age sample (12 of 13; 92.3%).

Considering the larger dataset of n = 98 individuals with presence/absence data for at least one stress indicator to maximize sample sizes, we observed a distinct difference between the Neolithic and Bronze Age patterns (Fig. 2). Specifically, porotic hyperostosis is common in the Neolithic sample (15 of 22; 68.2%), while linear enamel hypoplasia is relatively rare (7 of 25 individuals; 28.0%). The opposite is true in the Bronze Age, with 4 of 17 (23.5%) positive for porotic hyperostosis yet 13 of 17 (76.5%) with linear enamel hypoplasia.

We next tested whether the presence of paleopathological indicators of stress is predictive of individual-level deviations from the overall relationship between osteological stature and polygenic height score estimates. Based on the subset of individuals with presence/absence data for all three paleopathological indicators, the 41 individuals with one or more stress indicator were −0.917 cm shorter than expected on average compared with the 17 individuals without any stress indicators (). With the moderate magnitude of this difference and the relatively small sample size available for this analysis, this result was not unlikely based on chance expectations (t test; P = 0.555). The 18 individuals with two or more stress indicators were −1.175 cm shorter than expected on average relative to individuals with no stress indicators (P = 0.561).

Using our larger dataset, we next separately analyzed the effect of each paleopathological stress indicator on the osteological stature–polygenic height score relationship (Fig. 2). We found that the n = 21 individuals with cribra orbitalia were slightly but not significantly shorter (−1.23 cm) than expected on average compared with the n = 60 individuals without cribra orbitalia (P = 0.461; false discovery rate = 0.771). Linear enamel hypoplasia and porotic hyperostosis presence/absence were negligibly associated with osteological stature vs. polygenetic height score residual variation (Fig. 2 and ). These patterns did not change appreciably when excluding the few individuals with active cribra orbitalia or porotic hyperostosis lesions from their respective analyses (active lesions reflect adult stress while potentially masking evidence of lesions from childhood) ().

Finally, we investigated the relationship between cribra orbitalia, porotic hyperostosis, and linear enamel hypoplasia presence/absence and osteological vs. polygenic height score residual in the context of cultural period (). Of note, the n = 6 Neolithic individuals with cribra orbitalia were −4.12 cm shorter than expected on average compared with the n = 19 individuals from the same cultural period without cribra orbitalia (Fig. 2) (P = 0.329). This preliminary but suggestive effect was nearly absent in the Copper Age (+0.065 cm difference; P = 0.973), the only other period with sufficient presence/absence sample sizes for analysis.

Discussion

Bioarchaeologists have equated repeated observations of relatively shorter average adult statures in the Neolithic to a likely general health decline for individuals during this cultural period (13, 14, 16, 34, 55, 87). Combinations of reduced nutritional diversity, unpredictable food availability (e.g., crop failure, storage loss), and increased infectious disease burden may have negatively impacted childhood health and growth (18, 88, 89). Understandably, those prior studies did not account for the contribution of interindividual variation in the contribution of heritable genetic factors to adult stature. Yet, this consideration is especially important in light of updated understandings of considerable migration and gene flow processes associated with various farming transitions (90–93).

In our study, we sampled 167 prehistoric European individuals for whom both genome-wide ancient DNA data and intact long bones were available for analysis, making it possible to test whether Neolithic individuals were still osteologically shorter than expected when accounting (at least partly) for individual-level genetic contributions to height. Using this approach, we found that the average Neolithic farmer was indeed relatively shorter than expected compared with pre-Neolithic individuals (Fig. 1). Average osteological vs. expected stature then increased over each post-Neolithic cultural period. This gradual recovery may reflect a history of continued (although variable) cultural and technological innovations that ameliorated and/or overpowered the initial nutritional and disease stressors faced by the earliest farmers (64, 66–69, 94).

Our framework is related to but differs from that of a previous study by Cox et al. (58), who compared population-level osteological height estimates [n = 1,159 total individuals; the osteological data were from Niskanen et al. (95)] and ancient DNA-based polygenic height scores (n = 1,071) across prehistoric Europe. These two estimates were computed separately (i.e., typically not for the same individuals), thereby facilitating the large sample sizes. In contrast, our approach expressly considers individual-level dynamics in the relationship between these two variables, which is sample size restrictive yet potentially insightful. Interestingly, Cox and colleagues observed that mean osteological stature estimates and polygenic height scores were both similar between the European Mesolithic and Neolithic (58, 95). This result is in contrast to our own osteological height estimate observations in this study and those of previous bioarchaeological studies (15, 48), which again may reflect potentially interesting interpopulation variability as part of the nuanced complexity underlying subsistence shifts (94).

Limitations.

One important potential limitation of our study, aside from the number of individuals, is uncertainty regarding the ultimate portability of polygenic scores over genetic, geographic, and temporal distances (81, 82, 96, 97). Our analyses were also based on incomplete ancient DNA genotype data (however, note that we did exclude potential error from deamination-based ancient DNA damage by masking all potentially affected sites in the primary versions of our analyses). Yet, the significant, positive relationship between polygenic height scores and estimated osteological statures across our overall sample (Fig. 1) demonstrates the biological plausibility of our model. Moreover, our primary results were unchanged when we incorporated archaeological site latitude and longitude variables into our analyses.

However, even with complete genome-wide genotype data, polygenic height scores only capture a proportion of the heritable component of stature variation (83, 98–100). Furthermore, only a subset of variants is “available” for assessment in ancient genomic data. The predictive accuracy of polygenic scores is also higher when individuals are more closely related to the cohort used to develop the GWAS and generate the effect sizes of associated single nucleotide polymorphisms (SNPs), as increasing genetic distance from this cohort has been shown to negatively impact the predictive accuracy of polygenic scores (97, 100, 101). Therefore, for multiple reasons, our primary analytical approach might incompletely capture the stature-relevant effects of any genetic ancestry variation across our sample. That is, with respect to our hypotheses, we could only be partially accounting for any cultural period-confounded migration/gene flow among populations with different genetic height profiles. For example, gene flow as a result of the spreading of the Yamnaya/Corded Ware cultures (“steppe ancestries”) starting ∼4,000 to 5,000 y ago may have been associated with the introduction of relatively greater proportions of “tall” alleles into various regions of Europe (58, 64).

If this is a general phenomenon that extends to small effect and other loci not included in our individual-level polygenic height score calculations, then our cultural period–related inferences could be erroneous. Future improvements in the paleogenomics field related to low-coverage imputation and modeling the potential impacts of population turnovers and/or admixture events on polygenic scores will be important in driving forward genotype–phenotype trait reconstructions (86).

In the meantime, to help explore the potential effects of these processes on our results, we conducted a parallel set of analyses in which we included factors reflecting genome-wide genetic ancestries (from the first four components of a MDS analysis) in our model. When including both the polygenic height scores and MDS factors in the analysis, our downstream finding of shorter osteological statures than expected (i.e., based on polygenic height scores and MDS factors and in the context of the overall sample) in the Neolithic relative to other cultural periods was in fact attenuated. This approach does represent an overcorrection or more precisely, some level of double accounting of the genetic effect on height. That is, the MDS factors and the polygenic height scores are intertwined and cannot be deconvoluted readily due to the pervasive polygenicity of height (62, 102). Nonetheless, the signal attenuation result calls for cautious interpretations of the main results in our study while awaiting expanded future sample sizes.

To provide further insight into the relative magnitudes of these genetic components, we observed that 7.96% of the variation in the difference between individual and per-sex mean osteological stature is explained by individual-level polygenic height score alone (Fig. 1), whereas 8.4% of the osteological variation is explained by MDS factors alone (i.e., without polygenic height score included) (). Meanwhile, 14.2% of this variation is explained by a model including both individual polygenic height scores and the MDS components (). While these comparisons are imperfect for the reasons noted above, our polygenic height score does appear to contain valuable biological information for our hypothesis-testing approach.

Additionally, we note that comparable r2 values for the relationship between polygenic height score and estimated stature between those for individual modern human populations and our temporally and culturally diverse ancient sample are not necessarily expected (). That is, we have prior archaeological knowledge of major subsistence patterns and other cultural changes—with potentially large somatic effects on growth thereby reducing the strength of the relationship between polygenic height score and osteological adult stature—that occurred across our sample. How these cultural changes are associated with the (potentially highly variable) relationship between polygenic height scores and estimates of adult stature is the analytical essence of our study.

Toward the Integration of Paleopathological Data.

For a subset of the individuals in our study (n = 98), we were additionally able to consider the extent to which three paleopathological markers of nonspecific childhood and childhood-inclusive stress (linear enamel hypoplasia, cribra orbitalia, and porotic hyperostosis; each of which are maintained in the skeleton into adulthood) are associated with the relationship between osteological stature estimates and polygenic height scores. We observed a slight trend of relatively shorter than expected (given polygenic height score) adult statures among individuals with one or more childhood stress markers present ().

However, larger sample sizes will be necessary to more fully explore interplays among specific paleopathological indicators, osteological vs. genetic height scores, and cultural periods. Still, our findings do at least suggest that factors underlying skeletal growth trajectories are separable, at least in part, from those leading to paleopathological indicators of stress. In particular, high rates of paleopathology are still observed in post-Neolithic cultural periods (for example, in the Bronze Age), even after absolute osteological stature and actual vs. expected stature averages have recovered.

Looking further and more generally forward to future studies aiming to integrate per-individual osteological and ancient DNA information, we note that genetic risk for certain diseases that can lead to paleopathology may vary among individuals within a population and on average between groups. For example, Berens et al. (103) reported higher genetic risk scores for dental/periodontal disease in the agriculturalist relative to the pastoral populations they studied.

Summary.

We united previously disparate osteological and paleogenomic datasets for 167 prehistoric European individuals on a per-individual basis. Our results represent an advance in the study of whether and how a major cultural transition in human evolution affected physiological health. In particular, we show that the average Neolithic individual may have been relatively short even when correcting for expected individual genetic contributions to adult stature. This result may reflect reduced nutrition and/or increased infectious disease burden. We also preliminarily developed a framework for further consideration of these results in the context of particular paleopathological indicators of childhood stress. Looking forward, our model can be expanded in various dimensions (for example, to different world regions or to more constrained spatial and temporal contexts) in order to further the study of emergent physiological trade-offs across periods of dramatic cultural or environmental change. Integrated osteological–genetic approaches will increasingly become important components of the tool kit for studying the dynamics of past human health.

Materials and Methods

The 167 individuals in our dataset have broad geographical, temporal, and cultural period ranges. Radiocarbon or archaeologically calibrated dates, latitude/longitude coordinates, genetic sex, and archaeological/cultural period were obtained from the original paleogenomic and archaeological publications of the paleogenomic data used for this study (Dataset S1).

Processing Ancient DNA Sequence Reads.

Published sequence data were downloaded from the Sequence Read Archive and the European Nucleotide Archive databases as indicated per the respective papers alongside sequence data from in-process studies (Dataset S1). FASTQ (65, 67, 104–107) and BAM (Binary Alignment Map) (66, 68, 69, 108–122) files from shotgun sequence (n = 26) and DNA capture-based datasets (n = 141) were aligned/realigned using the Burrows–Wheeler aligner aln (123) to the human reference genome (hg19, build 37) with seeding disabled. For unaligned shotgun sequence data (n = 7) (65, 104–107), leeHom (124) was used to trim adapters and merge reads. Reads from sequence libraries that were not treated to remove damage signatures typical of ancient DNA (“non-uracil–DNA–glycosylase”) were subject to rescaling using mapDamage 2.0 (125) after mapping (105, 121), which downscales quality values for likely ancient DNA misincorporations based on read position and damage pattern. For partial and full uracil–DNA–glycosylase–treated libraries, two base pairs (bp) at the 5′ and 3′ ends were trimmed (prior to mapping) using seqtk (source code from ref.126) so as to not confound downstream analyses with potential postmortem deamination at the terminal read ends (66, 68, 109, 127). SAMtools was used to sort mapped reads and filter for mapping quality 30 and minimum bp 30, with duplicates removed using SAMtools rmdup (128). Read groups were added using Picard Tools AddOrReplaceReadGroups function (129). Following the GATK (Genome Analysis Toolkit) workflow (130), realigning indels was performed using RealignerTargetCreator and IndelRealigner, followed by BaseRecalibrator to minimize sequence error introduced by potential mismatches to the reference (131).

Genotyping and Imputation.

We implemented GATKUnifiedGenotyper (130) followed by imputation to maximize the amount of genetic information for downstream analyses and since genotype likelihood scores can be generated for imputation. We opted to impute diploid genotypes and missing sites for the individuals in our dataset [using the 1000 Genomes Project Consortium (132) reference panel] to minimize potential reference bias that may otherwise occur when using an alternative approach of randomly sampling one allele at each site (133).

The 1000 Genomes phase 3 genetic variants reference panel was used for genotyping and imputation, as provided by BEAGLE (134, 135). After removing variants that are not SNPs, multiallelic SNPs, and X and Y chromosomes, 77,818,345 variants remained. UnifiedGenotyper (130) was used to obtain genotypes and likelihood scores with the following parameters: –genotyping_mode –alleles <1000 Genomes reference panel>, –GENOTYPE_GIVEN_ALLELES, –output_mode EMIT_ALL_SITES, –AllSitesPLs -R (131). Due to the potential for postmortem damage impacting C > T and G > A allele changes, the per-chromosome VCF (Variant Call Format) files of the called genotypes were filtered for potential postmortem damage (modifying source code from ref. 136 for multiindividual VCF files and removing deamination in both directions). Potential deamination signal C > T (T > C) and G > A (A > G) genotypes were replaced in ancient individuals heterozygous for these genotypes with “./.” similar to previous paleogenomic studies (65, 76).

Genotype likelihoods were then estimated using the per-chromosome VCF files, followed by imputation of missing SNPs based on the genotype probability score using the 1000 Genomes phase 3 haplotypes (137) and GRCh37 genomic maps (138). Parameters for estimating genotype likelihoods were gprobs = true, gl ≤ input genotypes from UnifiedGenotyper> ref ≤ Beagle imputation reference panel>, map ≤ hg19 recombination map>. Imputation parameters were gt ≤ GL_output_VCF>, gprobs = true, impute = true, ref ≤ Beagle imputation reference panel>, map ≤ GRCh37 recombination map> (131). This resulted in 30,761,499 markers imputed/genotyped across 167 individuals. Prior to downstream analyses, the imputed VCF was filtered for a minimum genotype probability of 0.99 to maximize confident genotype calls postimputation and filter out less confident calls. We repeated the above pipeline of genotyping calling and imputation without filtering for potential deamination signals, and the results were consistent with the deamination filtered data ().

Assessing Accuracy of Imputation.

For the individuals in our study, the genomic/SNP capture sequence coverage mean for SNP positions analyzed was =1.595X, SD = 2.438, and median = 0.717, with a range of 0.002X to 22X. Prior to imputation, an average of 338,074 SNP sites per individual were covered by a minimum of one read with a median called genotype proportion =0.095 (i.e., median proportion of missing SNPs is =0.905) (). Our imputation process generally increased the number of called/inferred variants available for study (mean = 29,497,786; SD = 246,336; median = 29,591,650), reduced the proportion of missing SNPs (mean = 0.041; SD = 0.008; median = 0.038), and increased the proportion of called/inferred genotypes (mean = 0.959; SD = 0.008; median = 0.962) (). Specific to the genetic height scoring analysis, 1,658 SNPs were included, and of these markers, 155 were fully imputed variants (sites that were completely missing among all individuals).

A remaining limitation impacts paleogenomic studies in general, whether diploid genotypes were inferred via imputation or a “pseudohaploid” approach related to the broadly low coverage of ancient genomic data, which could impact the predictive accuracy of polygenic scores. To consider the potential effects of this phenomenon on imputation and genotype accuracy in our dataset, we evaluated how sequence coverage may have impacted our imputation accuracy (i.e., whether imputation is outperforming under high vs. low sequence coverage conditions). Specifically, we compared our imputed genotype data in the full dataset of the high-coverage Loschbour individual (∼16×) (115) with down-sampled BAM files (using SAMtools -s parameter) (128) from 3× coverage to 0.3× for chromosome 1 using SnpSift (139). We obtained a concordance rate (in terms of total sites recovered) of ∼97 to 99%, suggesting that imputation accuracy is not dramatically lower in the low-coverage imputed genotype data (). We also assessed the imputation accuracy of heterozygous sites by comparing the not imputed high-coverage genotype data for Loschbour with each down-sampled imputed BAM file from 0.3× to 3× coverage. At the lowest coverage of 0.3×, ∼85% of the heterozygous sites are still recovered with increase to ∼98% at 3× coverage ().

Polygenic Height Scores.

Polygenic scores were computed by downloading a publicly available GWAS dataset from the UK Biobank (77), specifically genome-wide summary statistics available from the Neale Lab (140). The quality controls implemented for the publicly available UK Biobank dataset (e.g., minor allele frequency > 0.1%, Hardy–Weinberg equilibrium P value > 1e-10 in 337,199 individuals) contain 10.8 million analyzable SNPs (141). The variant identifications (“variants.tsv.bgz”) were merged with the “standing height GWAS” file (“50_irnt.gwas.imputed_v3.both_sexes.tsv.bgz”), where the beta values represent the effect size of the “ALT” (alternate) allele subsequently used for performing the polygenic scores.

For our data, polygenic height scores were estimated using PLINK 1.9 (142) with clumping of independent SNPs. Clumping was used to identify the SNP with the lowest P value in each linkage disequilibrium (LD) block (142). This approach retains SNPs with the strongest statistical evidence while reducing the correlation between the remaining SNPs (143). Although all common SNPs could be used in polygenic scoring, clumping to remove SNPs that have limited statistical association is also a practical approach (144). Clumping was performed at the genome-wide P value 5e-08 using PLINK 1.9 parameters “–clump-r2 0.1” and “–clump-kb 1000” with the 1000 Genomes “Europeans” reference population panel to retain the most correlated SNPs (“index SNPs”) from the UK Biobank height summary statistics, which were then used to calculate the polygenic height scores. Polygenic scores were calculated using “–geno 0” to exclude missing genetic markers and the “–score” flag, extracting the specific index SNPs (131).

To help assess whether the level of portability between modern reference (UK Biobank) and ancient (“genetically distant”) individuals in our study was in line with expectations, we created approximately equivalent SNP and sample size datasets for individuals who were excluded from the UK Biobank cohort (“non-GWAS” individuals; conducted using the UK Biobank Resource under application no. 31063 to R.M.). These excluded individuals have differing genetic relatedness from the UK Biobank cohort considered in the stature GWAS, the data that served as the basis for our per-individual polygenic height score estimates. This non-GWAS UK Biobank cohort is a very rough but conceptually similar approximation to the sample of ancient individuals analyzed in our study.

Specifically, we evaluated the predictive accuracy between phenotypic and predicted genetic height (r2) for two sample sets (n = 361,182 GWAS and n = 4,712 non-GWAS individuals), each stratified by sex and for both sexes combined. A dataset of 784,256 genotyped SNPs was subject to LD clumping with parameters similar to those in our ancient DNA analysis (using PLINK 1.9). Postclumping, 5,183 SNPs were retained at the genome-wide significance level. We also performed a down sampling (100 random draws) of the full set of GWAS and non-GWAS cohorts to 167 individuals.

For females and males combined, the non-GWAS individuals exhibit slightly lower variance in the relationship between phenotypic height and predicted genetic height (r2 = 0.031) relative to the GWAS cohort (r2 = 0.04) (). For the down-sampled data (167 individuals), there is similarly a slightly lower variance in the relationship between stature differences and polygenic height scores among non-GWAS individuals (r2 = 0.037) relative to the GWAS cohort (r2 = 0.043). Relative to these results, the variance explained by polygenic height score in our ancient sample (r2 = 0.0796) is not exceedingly low.

Long Bone Measurement and Stature Reconstruction.

Both newly collected (n = 93 individuals) and previously collected/published (n = 54) osteological data were included in this study (Dataset S1). For n = 20 of the latter set of individuals, only precalculated terminal height estimates (based on unavailable long bone length measurements) were available (Dataset S1). Only adult individuals were included in our study. For newly collected data, this assessment was based on the complete fusion of all long bone epiphyses; we otherwise relied on classifications of “adult” in the published record (which were likely based on the same criterion).

Permissions to collect new long bone measurement data were coordinated with researchers (coauthors on this publication) affiliated with the museums and university departments housing the various individuals. An osteometric board was used to measure the maximum length measurements (to the nearest millimeter) of the femur, tibia, radius, and humerus following standard osteological methods (44, 144, 145). Intact long bones were selected, either the left or right side depending on availability and preservation; if both sides were available and fully preserved, then both were measured.

We used a regression-based approach to reconstruct osteological stature from Ruff et al. (44). These equations were developed using 501 individuals from across Europe ∼7,000 BC to AD 1900, broadly approximating the geographical and temporal span of the individuals in our dataset. Sex-specific regression equations for the femur, tibia, humerus, and radius were used (), with SE estimates ranging from 1.66 to 2.73% (44). For the tibia, separate “north” and “south” equations are available (44). Given the potential for migration occurring across the temporal and cultural periods in our dataset, we computed estimates from both equations and averaged them for the tibia-derived stature estimate for all individuals in our study (with available tibia measurements) regardless of geographic origin. When measurements from multiple different bones from the same individual were available, stature estimates derived from each of the different bones were estimated separately and then averaged to obtain a single point estimate per individual.

For paleopathological evaluations, 82 individuals were newly characterized, and 16 were published/previously characterized (). Crania with at least one permanent incisor were examined to record the presence or absence and severity of three skeletal indicators of nonspecific stress: porotic hyperostosis, cribra orbitalia, and linear enamel hypoplasia (21, 22, 27, 28). Cribra orbitalia was assessed according to Stuart-Macadam (21) on a five-stage scale of severity (n = 81 evaluated) and whether the lesions were healed or active (n = 80) (21). Porotic hyperostosis was evaluated on a three-stage scale (n = 80) and whether the lesions were healed or active (n = 79) (22), and linear enamel hypoplasia was assessed as present or absent (i.e., whether one or more linear bands of decreased enamel thickness were visible; n = 81) (27). The differences in the number of individuals assessed for healed or active lesions are because of the exclusion of a previously published individual (n = 1 Neolithic) (29) identified as having active lesions and n = 1 Mesolithic individual for whom the nature of the porotic hyperostosis lesions was unspecified (146).

Statistical Analyses.

Statistical analyses were conducted in RStudio (v1.2.5033). Our main linear model was generated using osteological height and genetic height scores with sex as a covariate. Data normality was assessed using “ggResidpanel” (v0.3.0) (147) (). The residuals from this model were the basis for downstream analyses comparing patterns of stature variation across cultural periods as well as with the paleopathology data. Statistical analyses (t tests) on the residuals from various linear models (including those below) were performed, and the results are provided in .

We performed two additional analytical iterations of our linear model to evaluate consistency of downstream results. First, we included latitude and longitude as additional factors in the linear model framework described above. Second, we included factors related to genetic ancestries variation into our linear model following the approach of Cox et al. (84) using MDS rather than PCA.

Following Cox et al. (84), we used an MDS-based rather than a PCA-based approach to account for population stratification. PCA uses genetic correlations among individuals to represent population structure (148, 149), while MDS converts population structure into a matrix of observed (pairwise) genetic distance among individuals (based on similarity/dissimilarity) (85, 142, 150). Essentially, MDS configures the coordinates of the data in a low-dimensional space in a way that reflects the pairwise relationships of the original data, which is appropriate for large genomic datasets, particularly when the data are nonlinear (i.e., missing information) (151). In this way, MDS (152) reduces the “noise” from the combined influences of genetic variation and phenotypic association on height by accounting for the potential risk of false positives due to residual population stratification from individual ancestries (85).

Specifically, we performed the MDS analysis in PLINK (v1.9) (142) using the full genome-wide SNP genotype data available for all of the individuals in our study. We generated the “plink.genome” file (plink –file –genome), which was used as input for the MDS analysis specifying four dimensions (plink –file –read-genome plink.genome –cluster –mds-plot 4). The first four components (C1, C2, C3, and C4) were each included in the linear model.

Since the selection of the number of MDS components is somewhat arbitrary (but following the methods of another study) (84), we also fit increasing numbers of components (from five to eight) and then performed all subsequent analyses accordingly. In each case, the signal of decreased but muted patterns of reduced stature among Neolithic individuals relative to other cultural periods is replicated (e.g., Neolithic average residual for MDS 8 = −0.589 ± 6.8 cm vs. −0.603 ± 6.9 cm for MDS 5) ().