Polygenic Selection and Environmental Influence on Adult Body Height: Genetic and Living Standard Contributions Across Diverse Populations

Abstract

We analyzed whole-genome sequencing (WGS) data from 51 populations and combined WGS and array data from 89 populations. Multiple types of polygenic scores (PGS) were employed, derived from multi-ancestry, between-family genome-wide association study (GWAS; MIX-Height), European-ancestry, between-family GWAS (EUR-Height), and European-ancestry siblings GWAS (SIB-Height). Our findings demonstrate that both genetic and environmental factors significantly influence adult body height between populations. Models that included both genetic and environmental predictors best explained population differences in adult body height, with the MIX-Height PGS and environmental factors (Human Development Index [HDI] + per capita caloric intake) achieving an R² of .83. Our findings shed light on Deaton’s ‘African paradox’, which noted the relatively tall stature of African populations despite poor nutrition and childhood health. Contrary to Deaton’s hypotheses, we demonstrate that both genetic differences and environmental factors significantly influence body height in countries with high infant mortality rates. This suggests that the observed tall stature in African populations can be attributed, in part, to a high genetic predisposition for body height. Furthermore, tests of divergent selection based on the Q_ST (i.e., standardized measure of the genetic differentiation of a quantitative trait among populations) and F_ST (neutral marker loci) measures exceeded neutral expectations, reaching statistical significance (p < .01) with the MIX-Height PGS but not with the SIB-Height PGS. This result indicates potential selective pressures on body height-related genetic variants across populations.

Human body height (stature) is a complex trait that exemplifies the intricate interplay between genetic predispositions and environmental influences. Twin and family studies have consistently demonstrated the high heritability of body height, with estimates ranging from 80−90% in Western adult populations (Grasgruber et al., 2014; Silventoinen et al., 2003). However, this strong genetic component presents a paradox when juxtaposed against the marked increase in average body height observed globally over the past century, largely attributed to improvements in living conditions.

Heritability, in this context, refers to the proportion of variance in a trait within a population that can be attributed to genetic differences. High heritability means that genetic differences play a substantial role in determining individual differences in body height. However, it is crucial to understand that even a heritability of 100% does not preclude environmental influences on a trait in the more general sense. For instance, a new kind of environmental cause may be introduced to change the phenotype. If the distribution of this cause varies between people, the heritability will be decreased (total variance increases while genetic variance remains constant). Or, it may leave heritability unchanged if the causal factor is evenly distributed and lacks interactive effects with genetic causes (total variance remains the same, but the mean is changed). Because of this, heritability estimates can vary across different environments and populations; they are not strictly speaking a property of the trait itself, but a population statistic that depends on the context. The large increase in human height seen over the last few hundred years, then, is not a paradox because it shows that while genetic causation dominates individual differences within a cohort of a population, environmental factors between cohorts may have a large effect.

While we do not know for certain what the causes of the historical increase in height are, nutrition, healthcare access and overall socioeconomic development have been linked (NCD Risk Factor Collaboration, 2020). Specific factors include protein and overall energy (caloric) intake, disease prevalence, infant health, and living conditions such as access to clean water and sanitation (Checkley et al., 2008; Martorell & Zongrone, 2012; Prüss-Üstün et al., 2014; Victora et al., 2008).

These environmental factors influence body height through various biological mechanisms, as must be the case for any variable affecting humans. Adequate nutrition, particularly protein intake, provides the building blocks necessary for bone growth and development. Caloric intake ensures sufficient energy for growth processes. Disease prevalence in childhood can impact body height by diverting energy from growth to immune responses and by interfering with nutrient absorption. Early life health conditions, including maternal health during pregnancy and infant nutrition, set the stage for future growth trajectories. As such, access to clean water and sanitation reduces the risk of infections that could impair growth. Together, these factors create an environmental context that either supports or hinders an individual’s potential for body height.

Height’s sensitivity to environmental conditions has made it a valuable proxy for living standards in economic and social science research. Variations in body height across populations and time periods often reflect cumulative impacts of nutrition, healthcare and socioeconomic conditions during growth periods, with taller stature generally associated with better childhood conditions and improved later-life outcomes (Deaton, 2007). Using body height as an indicator of development between populations or across long time spans, however, does rely on a blank slatist assumption that populations do not differ in their genetic potential for body height, and that these do not change over time (no selection or genetic drift). This assumption is in question in the light of persistent height gaps between modern populations, even when these grow up in the same countries, and because studies of both modern and ancient genomes show that height was under selection (Piffer & Kirkegaard, 2024; Stulp et al., 2015; Stulp et al., 2023).

To investigate these dynamics, we employed regression models to assess the impact of height polygenic scores (PGSs) and environmental factors on body height across different countries and ethnic groups. In the field of evolutionary ecology, the potential association between populations’ phenotype and their local environmental conditions is a crucial step in identifying the selective pressures responsible for adaptive phenotypic differentiation (Blanco-Sánchez et al., 2024; Blanquart et al., 2013). We utilized the Human Development Index (HDI) as a general indicator of living standards, infant mortality rates to indicate infant health, and measures of daily per capita protein and total energy intake, along with the prevalence of wasting, as indicators of nutrition.

PGSs are numerical estimates of an individual’s genetic predisposition or potential for a trait, based on the total effect of many genetic variants identified through genome-wide association studies (GWAS). These scores are typically calculated by summing the number of trait-associated alleles an individual carries, weighted by the effect size of each allele. PGSs provide a way to quantify the genetic component of complex traits like body height, allowing for the investigation of genetic influences across populations.

This study sought to assess the validity of Deaton’s (2007) hypothesis, which suggests that the positive correlation between health, nutrition and average height observed in affluent nations does not hold true for developing countries. Furthermore, we investigated the paradoxical phenomenon highlighted by Deaton, wherein African populations exhibit relatively tall statures despite poorer childhood health and nutritional conditions. Our research aimed to determine whether genetic factors could account for this unexpected trend.

Genetic variation across populations often exhibits spatial patterns, with geographically proximate groups showing greater genetic similarity. This phenomenon, known as spatial or genetic autocorrelation, can be quantified using measures such as the fixation index (F_ST), which assesses population differentiation (genetic distance) based on allele frequencies. By incorporating these metrics, we can account for shared ancestry and migration patterns, enabling a more nuanced analysis of the relationship between genetic factors and body height across diverse populations. This approach also facilitates the identification of signals indicating local adaptation (Berg & Coop, 2014). To further investigate potential selective pressures, we employed two standard population genetics tests for divergent selection: the Q_ST-F_ST comparison (Spitze, 1993) and the F_ST-enrichment test (Guo et al., 2018). Significant differences between population genetic differentiation in phenotypic traits (Q_ST) and population differentiation in molecular markers driven only by neutral processes (F_ST) are interpreted as evidence of selection (Leinonen et al., 2013; Merilä & Crnokrak, 2001; Whitlock, 2008; Whitlock & Guillaume, 2009.

By integrating diverse genetic datasets with detailed environmental data, this study aims to provide a comprehensive understanding of how genetic predispositions and environmental contexts interact to determine body height variations globally. This research seeks to contribute to the broader discourse on human growth and development by elucidating the relative contributions of genetic and environmental factors.

Methods

Genome Samples

Samples were collected from a variety of public databases and were classified by the technology employed to genotype the DNA (micro-array chips vs whole-genome sequencing (WGS).

The WGS dataset was comprised of the following: 1000 Genomes (1000 Genomes Project Consortium, 2015), gnomAD (Karczewski et al., 2020), SweGen (Rentoft et al., 2019), Genome of the Netherlands (GoNL, nlgenome.nl), 1000 Polish genomes (Kaja et al., 2022), NARD (Yoo et al., 2019), Turkish Genome Project (Alkan et al., 2014), Taiwan Genomes (Hsu et al., 2023), Mexico City Prospective Study (Ziyatdinov et al., 2023), Korean Variant Archive (Lee et al., 2022), ABraOM (Naslavsky et al., 2022), DanMAC5 (https://danmac5.cpr.ku.dk/) covering 51 ethnic groups in total.

The array dataset was compiled from a variety of open access online resources. We imputed the datasets after QC using TOPMed imputation server (Das et al., 2016) with Minimac imputation (Fuchsberger et al., 2014), covering 89 ethnic groups in total.

GWASs Used for Polygenic Scores

We utilized results from multiple recent GWASs in order to test the robustness of the results. Each model was named after the population it was trained or, or the method used.

1. 1. MIX-Height. For body height, we used the significant single nucleotide polymorphisms (SNPs) from the largest GWAS to date (Yengo et al., 2022), which comprised a multi-ancestry sample (after LD pruning with a threshold of r ² < .1).

2. 2. EUR-Height. We used the significant SNPs from the European-ancestry subsample of the Okbay et al. (2022) GWAS.

3. 3. SIB-Height. Summary statistics for sibship (within-family) GWAS of body height were retrieved from the largest meta-analysis of sibship GWAS (Howe et al., 2022). There were 290 SNPs that remained significant after clumping and applying a significance threshold (p < 5 × 10-8, LD r ² < .1). Within-family GWASs, specifically sibship studies, offer unique advantages in genetic research. By comparing siblings who share approximately 50% of their genetic material and a common family environment, these studies can better isolate the effects of specific genetic variants on body height. This approach helps control for confounding factors such as population stratification and shared environmental influences that can bias traditional GWAS results. However, sibship studies also have limitations, including smaller sample sizes and reduced statistical power due to the lower phenotypic and genotypic variance within families compared to between families. Despite these practical constraints, sibship GWAS provide valuable insights into the genetic architecture of complex traits like body height (Howe et al, 2022).

Health, Nutrition and Living Standards

Food and Agriculture Organization (FAO) daily per capita protein supply estimates were obtained from Our World in Data (https://ourworldindata.org/grapher/daily-per-capita-protein-supply) and the average infant mortality (from birth to 1-year-olds only) prevalence by country for the years 1995−2000 was obtained from the United Nations 2022 Revision of World Population Prospect (https://population.un.org/wpp/).

Wasting prevalence (the share of child under 5 years old that fall two standard deviations below the expected weight for their height) was obtained from the estimates provided by the World Bank (https://ourworldindata.org/wasting-definition) for the years 1995−2000.

For HDI, we used the estimates provided by the United Nations Development Programme (2024). We obtained subnational HDI estimates for Italy, Spain, China and the US from www.globaldatalab.org.

Average daily dietary energy consumption per capita for 2018 were obtained from the report published by the FAO (Roser et al., 2013).

Phenotypes

Average measured body height (the phenotype, or trait) was obtained mainly from the largest meta-analysis by the NCD-RisC factor collaboration (2020). Other sources were used for subregions that were not covered by the latter. The body height estimates for the Chinese provinces were obtained from Piffer and Kirkegaard (2024), who relied on two sources: (1) an analysis of 57,574 samples from the Chinese General Social Survey (CGSS; Lu et al., 2022), and (2) an analysis based on data from 660K users collected from the big data platform ‘Xiangshan Weighing Instrument Group’ (2022年中国居民身高体重健康数据报告 [2022 Height and Weight Health Data Report of Chinese Residents]). For Sudan Nilotic, Chali (1995) was used. For Italy, we obtained regional estimates from Corsini (2009).

Spatial Autocorrelation

Spatial autocorrelation measures how a variable correlates with itself across geographic locations, determining whether nearby observations are more similar (positive autocorrelation) or dissimilar (negative autocorrelation) than expected by chance.

A typical example of positive spatial autocorrelation is housing prices, where high-value properties cluster together, as do lower priced ones, driven by shared factors like neighborhood amenities and socioeconomic conditions. Negative spatial autocorrelation can be seen in the placement of competing retail stores, which tend to spread out to avoid market saturation, placing distance between themselves and direct competitors.

In population genetics, spatial autocorrelation is often assessed by comparing geographic distances with allele frequencies or genetic distance measures like F_ST (Sokal et al., 1989). Piffer (2015) adapted this methodology, using F_ST distances as a proxy for spatial similarity, replacing geographic distances, and focusing on the absolute differences in PGSs to assess similarity at loci under selection. This innovative approach allowed for the first-time testing of autocorrelation purely at the genetic level, and by combining it with the partial Mantel test, it quantified how well PGSs predicted phenotypic differences while controlling for neutral genetic variation.

Mantel Test

To assess the relationship between genetic and phenotypic distances across populations, we employed the Mantel test (Hubert et al., 1981; Mantel, 1967; Sokal, 1979). This nonparametric method evaluates the correlation between two distance matrices.

The Mantel test calculates a z statistic, defined as:

$${Z_{AB}} = \sum\nolimits_{ij} {{a_{ij}}{b_{ij}}} $$

Where:

X_ij and Y_ij are elements of the two distance matrices X and Y, respectively.

This z statistic is then compared to its distribution under the null hypothesis of no correlation between the matrices. The significance of z is assessed through a permutation procedure:

The observed z statistic is calculated for the original matrices.

One matrix is randomly permuted, and z is recalculated.

This process is repeated numerous times to generate a null distribution.

The p value is determined as the proportion of permuted z statistics that are equal to or more extreme than the observed z statistic (Mantel, 1967; Mielke, 1979).

We used this test to evaluate the correlation between genetic distances (as measured by F_ST) and differences in body height or PGSs between populations. A significant positive correlation would indicate that genetically distant populations tend to have more divergent body height or PGSs, potentially suggesting a role for population structure or local adaptation in body height variation.

Moran’s I

To assess the spatial autocorrelation of body height and PGSs across populations, we employed Moran’s I statistic (Moran, 1950). Moran’s I is a measure of global spatial autocorrelation that quantifies the degree to which a variable is correlated with itself through space. The statistic ranges from −1 to +1, where:

Values of +1 indicate perfect positive spatial autocorrelation (similar values cluster together).

Values of −1 indicate perfect negative spatial autocorrelation (dissimilar values cluster together). Values near 0 suggest random spatial distribution.

Moran’s I is calculated as:

$$I = {N \over W} \cdot {{\sum\nolimits_i {\sum\nolimits_j {{w_{ij}}({x_i} - \overline x)({x_j} - \overline x)} } } \over {\sum\nolimits_i {{{({x_i} - \overline x)}^2}} }}$$

Where:

N is the number of spatial units; x is the variable of interest; ${\rm{\overline x}}$ is the mean of x; w_ij is the spatial weight between locations i and j; W is the sum of all spatial weights.

We used a distance-based weight matrix, where weights were inversely proportional to the geographic distance between populations. Statistical significance was assessed through a permutation test with 999 randomizations. This analysis allows us to determine whether the spatial distribution of body height and PGSs exhibits significant clustering or dispersion across the studied populations.

Spatial Autoregressive and Spatial Error Model

The spatial autoregressive model (SAR) incorporates spatial dependence directly into the dependent variable. It assumes that the value of the dependent variable in a given location is influenced by the values of the dependent variable in neighboring locations. This spatial dependence is captured through a spatial lag term, which represents the average values of the dependent variable in the neighboring regions, weighted by the spatial proximity of those regions.

The SAR model is formulated as follows:

$${\rm{Y = \rho WY + X\beta + \varepsilon Y}}$$

Where:

• Y is the dependent variable.

• ρ is the spatial autoregressive coefficient, which captures the strength of the spatial dependence.

• W is the spatial weights matrix, which defines the structure of spatial relationships (e.g., neighbors based on geographic or genetic distance).

• Xβ represents the explanatory variables and their coefficients.

• ε is the error term.

In this model, ρ indicates the degree to which the dependent variable in one location is influenced by the dependent variables in neighboring locations. The model can be used to measure both direct and indirect (spillover) effects of the explanatory variables. The SAR model is useful when there is reason to believe that the outcome at one location is influenced by outcomes in nearby locations. The spatial error model (SEM) accounts for spatial dependence in the error terms, rather than in the dependent variable itself. This is based on the idea that unobserved factors, which influence the dependent variable, are spatially correlated. In this case, spatial autocorrelation arises from omitted variables or measurement errors that are spatially clustered, rather than the dependent variable being directly influenced by neighboring values.

The SEM is formulated as follows:

$$Y = X\beta + \epsilon $$

With the error term specified as:

$$\epsilon = \lambda W\epsilon + u$$

Where:

• Y is the dependent variable.

• Xβ represents the explanatory variables and their coefficients.

• ε is the error term, which exhibits spatial autocorrelation.

• λ is the spatial autoregressive parameter for the errors, capturing the extent to which spatially correlated unobserved factors influence the outcome.

• W is the spatial weights matrix, similar to the SAR model.

• u is a white noise error term.

The SEM is useful when there is spatial autocorrelation in the residuals, indicating that unobserved spatial processes are influencing the outcome. This model helps correct for this spatial autocorrelation, improving the overall model fit and accuracy.

The following steps were employed to implement the spatial autoregressive and spatial error models:

1. 1. Defining F _ST thresholds. The thresholds variable defined different F_ST values (e.g., 0.03, 0.05, 0.07, 0.09). These thresholds determined which populations were considered neighbors based on genetic distance.

2. 2. Filtering neighbors by F _ST . The genetic distance matrix was filtered for each threshold to retain only neighbors with F_ST values less than the specified threshold. This ensured that only populations with sufficient genetic similarity (i.e., below the threshold) were treated as neighbors.

3. 3. Neighbor list creation. After filtering the distance matrix, a neighbor list was created for each population, including only those populations that met the F_ST threshold.

4. 4. Spatial weights list. The neighbor list was converted into a spatial weights list, representing the influence of each population on its neighbors. This spatial weights list was used as input to the spatial models.

5. 5. Model and impact calculation. For each F_ST threshold, a spatial autoregressive model was fitted using the spatial weights list. The direct, indirect, and total impacts of the spatial lag were then calculated, and the results were analyzed.

This process was repeated for each F_ST threshold to assess the influence of genetic similarity on the spatial relationships between populations.

Tests of Divergent Selection

Q_ST was computed using the formula Q_ST = σ²B/(σ²B + 2σ²W) (Leinonen et al., 2013). Q_ST is defined as the level of genetically based population differentiation in quantitative traits (Li et al., 2019).

The total (additive) genetic variance is the variance of the PGSs across all individuals in all populations. The genetic variance within populations is the average variance of the PGSs within each population, weighted by the number of individuals in each population.

Q_ST is then calculated as the genetic variance among populations divided by the sum of the genetic variance among populations and twice the genetic variance within populations.

In fact, population differences in the mean of a quantitative trait due to positive covariances — that is, (cross-population) linkage disequilibrium (LD) — between distant variants can arise for polygenic traits under divergent selection. In practice, this happens when alleles with similar effects are driven to similar frequencies within populations across multiple loci (Latta, 1998; Le Corre & Kremer, 2003; Ma et al., 2010). The other component of genetic differentiation in quantitative traits, F_ST, does not take into account this covariance of allelic effects and it was shown to be usually very small for highly polygenic traits subject to divergent selection pressures (Berg & Coop, 2014).

Indeed, F_ST is based on the variances of individual allele frequencies, which are unsigned, meaning that positive and negative effects of alleles on a trait can cancel each other out, leading to an underestimation of the true extent of genetic differentiation in polygenic traits.

Conversely, there may be substantial levels of genetic differentiation (F_ST > 0.15) without any variations in the population means (Q_ST = 0; Le Corre & Kremer, 2012).

If Q_ST/F_ST quantity is higher than 1, and p < .05, the quantitative trait in question is inferred to have been subject to divergent selection.

To produce Q_ST values free of the cross-population (long-range) LD, the effect and noneffect GWAS alleles were randomly shuffled with a probability of .5 to produce a null distribution of PGSs and calculate random Q_ST values. The expected value of Q_ST is equal to the F_ST of GWAS SNPs (F_STQ; Le Corre & Kremer, 2012) because this operation removes the variance due to cross-population linkage disequilibrium. We filtered variants with MAF < 0.01 (minor allele frequency) in any of the 5 1KG superpopulations. By filtering out variants with MAF < 0.01, the analysis focuses on variants that are more reliably measured, better understood, and more consistent across populations, thereby producing more robust and interpretable Q_ST estimates that align with the theoretical expectations set forth by studies like Le Corre and Kremer (2012).

Results

Tables 1 and 2 present the following data for each GWAS: the number of independent SNPs after linkage disequilibrium (LD) clumping and thresholding (r² < .1 and p < 5 × 10⁻⁸); the number of matching SNP IDs found in the WGS and WGS + array datasets; the Cronbach’s alpha of the resulting PGSs.

Table 1. WGS dataset

Note: WGS, whole-genome sequencing; GWAS, genomewide association study; SNP, single-nucleotide polymorphism.

Table 2. WGS + array dataset

Note: WGS, whole-genome sequencing; GWAS, genomewide association study; SNP, single-nucleotide polymorphism.

The PGSs computed using the WGS dataset consistently demonstrated higher match rates and internal consistency compared to those derived from the WGS + array dataset. This difference in performance is likely due to the more comprehensive genomic coverage provided by WGS.

Among the different GWAS sources, the multi-ancestry GWAS (MIX-Height) yielded PGSs with the highest internal consistency, as measured by Cronbach’s alpha. Specifically, the MIX-Height PGS achieved alpha values of .85 and .49 in the WGS and WGS + array datasets respectively. In contrast, the sibship GWAS (SIB-Height) produced PGS with lower internal consistency, with alpha values of .51 and .21 in the WGS and WGS + array datasets respectively.

These results suggest that the choice of GWAS source and genomic data type (WGS vs. WGS + array) significantly impacts the reliability and consistency of the resulting PGSs. The superior performance of the multi-ancestry GWAS-derived scores may indicate its greater generalizability across diverse populations.

Correlation of Polygenic Scores With Phenotypic and Environmental Variables

Our analysis revealed consistent patterns of correlation between average PGSs, environmental variables, and phenotypic body height across both the WGS and WGS + array datasets (Figures 1 and 2).

WGS Dataset

Environmental factors positively associated with body height included daily per capita energy intake, Human Development Index (HDI), and daily per capita protein intake, with correlations of 0.73, 0.70, and 0.62, respectively. Conversely, wasting prevalence and infant mortality were negatively correlated with body height, with correlations of −.61 and −.59 respectively. PGSs also showed positive correlations with average body height, ranging from .26 for SIB-Height to .70 for EUR-Height-2014.

WGS + Array Dataset

Environmental factors positively correlated with body height included daily per capita energy intake, HDI, and daily per capita protein intake, with correlations of .60, .56 and .50 respectively. Wasting prevalence and infant mortality were negatively correlated with body height, with correlations of −.52 and −.53 respectively. PGSs were also positively correlated with average body height, ranging from .34 for EUR-Height to .70 for EUR-Height-2014.

These results highlight consistent patterns across both datasets, with environmental factors and PGSs showing strong correlations with average body height.

Figures 1a and 1b (WGS and WGS + array dataset respectively). Correlation matrix of PGS means and phenotypic means.

Figure 1a. Whole genome sequencing (WGS) dataset: correlation matrix of polygenic score means and phenotypic means.

Note: *p < .05, **p < .01, ***p < .001.

Figure 1b. Whole genome sequencing (WGS) + array dataset: correlation matrix of polygenic score means and phenotypic means.

Note: *p < .05, **p < .01, ***p < .001.

Figures 2a and 2b display scatterplots visualizing the correlations between MIX-Height and SIB-Height with average body height within the WGS and WGS + array datasets respectively.

Figure 2a. Whole genome sequencing (WGS) dataset with scatterplots of MIX-Height, SIB-Height and average measured body height.

Figure 2b. Whole genome sequencing (WGS) and array dataset with scatterplots of MIX-Height, SIB-Height and average measured body height.

Regression Models

Table 3 summarizes the results of 30 models examining the relationship between PGSs, environmental factors, and body height across the two datasets (WGS and WGS + array). Per capita (log) GDP was excluded due to its high correlation (.95−.96) with HDI, which is a superior indicator of overall socioeconomic development.

Table 3. Regression models

Note: WGS, whole-genome sequencing; PGS, polygenic score; HDI, Human Development Index;

***p = .001 **p = .01 *p = .05.

Key findings include:

1. 1. Models using MIX-Height PGS generally show strong associations with body height, with PGS beta values ranging from 0.48 to 0.60.

2. 2. Environmental factors like HDI, calories, and proteins also show strong positive associations with body height, with beta values typically above 0.5 and often significant (indicated by ***).

3. 3. Combined models (e.g., HDI and calories) often show higher model R² values, indicating better explanatory power.

4. 4. Models using SIB-Height PGS show similar patterns but generally have lower PGS beta values and model R² values compared to MIX-Height PGS models.

5. 5. The WGS dataset generally yields higher model R² values compared to the WGS + array dataset.

Overall, the MIX-Height PGS and HDI, calories, and proteins as environmental factors are consistently significant predictors of body height, with combined models providing the highest explanatory power.

Effects of Environmental and Genetic Factors Among High Infant Mortality Countries

We subset the WGS + array dataset to only those countries or regions with high infant mortality rates, resulting in a subsample of 27 populations. Our regression analysis incorporated both genetic and environmental predictors, specifically PGSs for body height, infant mortality, and caloric intake.

The results presented in Table 4 demonstrate that both MIX-Height and SIB-Height PGS significantly predict body height, with MIX-Height showing stronger associations overall. Models that included infant mortality as an environmental predictor indicated a significant negative relationship with body height, while those with caloric intake showed a significant negative association as well.

Table 4. Regression models (only high infant mortality countries)

Note: WGS, whole-genome sequencing; PGS, polygenic score.

***p = .001 **p = .01 *p = .05.

Interaction Effects

We coded HDI as a dummy variable, with a cut-off value of 0.8 separating low from high HDI populations. The interaction between MIX-Height and HDI was not significant, but the one between Calories and HDI was significantly negative, indicating, contrary to Deaton’s predictions, that the positive effect of calorie intake on body height is reduced in the higher HDI group (Table 5).

Table 5. Regression models (with interaction term)

Note: WGS, whole-genome sequencing; PGS, polygenic score.

However, in the model with SIB-Height, neither the interaction between SIB-Height and HDI nor the one between Calories and HDI were significant (Table 5).

Latitude

Latitude showed a tendency towards positive correlations with body height PGS in both the WGS and WGS+Array datasets. However, statistical significance was limited. In the WGS dataset, only EUR-Height reached significance. These observations, combined with Bergmann’s rule linking body size to latitude (Bergmann, 1848), suggested the need for further investigation into latitude’s influence on both phenotypic and genotypic body height.

To explore patterns of local adaptation, we developed regression models using body height PGS as the dependent variable and latitude as the independent variable. We enhanced the model by including superpopulation as a categorical variable.

Although human genetic variation can be clustered in many ways, we defined six superpopulations, expanding on the five categories used in the 1000 Genomes Project (1KG): AFR (African); AMR (Amerindian); EAS (East Asian); EUR (European); SAS (South Asian); MENA (Middle Eastern/North African) — added to cover populations in our dataset not represented in 1KG.

Notable characteristics of our categorization:

• The AMR category encompassed both indigenous populations (e.g., Amazonian, Andean, Mesoamerican) and admixed groups (e.g., Mexican, Puerto Rican).

• Similar to 1KG, our AFR superpopulation included admixed samples, such as African Caribbean (ACB) and African ancestry in Southwest USA (ASW).

• Some heavily admixed populations could not be classified under this scheme; for example, Basters or Coloured from South Africa, or Brazilians.

This approach allowed us to account for both geographical and genetic diversity in our analysis of height-related genetic adaptations across different populations.

This allowed us to assess the impact of latitude while accounting for broad ancestry groups.

Controlling for superpopulation, latitude had a significantly positive relationship (beta from 0.31 to 0.52) with body height PGS in 3 out of 4 models. The effect of superpopulation, especially African and East Asian ancestry, suggests that the correlation between latitude and body height PGSs is obscured by ancestry (Table 6).

Table 6. Regression of body height polygenic scores on latitude and superpopulation

Note: PGS, polygenic score; AFR, African; EAS, East Asian; EUR, European; MENA, Middle Eastern/North African; SAS, South Asian; AMR (not shown because used as reference group), Amerindian/Hispanic.

***p = .001 **p = .01 *p = .05.

Spatial-Genetic Autocorrelation Models

The spatial autocorrelation tests were conducted using the combined WGS and array dataset. This dataset was chosen because it allowed for the calculation of F_ST (Fixation Index) genetic distances across all samples. In contrast, the WGS dataset alone contained numerous samples with frequencies derived from aggregate data, making F_ST calculations unfeasible for those particular samples.

Typically, spatial autocorrelation in the dependent variable (or residuals) is used to check whether the model adequately captures spatial dependencies in the outcome of interest. However, since we are also interested in the relationship between neutral genetic variation (as captured by F_ST) and variation at quantitative trait loci (i.e. PGS), we also computed the autocorrelation at the level of the PGS.

F_ST was computed using PLINK 2.0 (www.cog-genomics.org/plink/2.0/; Chang et al., 2015), with the default (Hudson) method (Hudson et al., 1992). After filtering for missing genotype (mind > 0.5), 251 samples removed due to missing genotype data, including two entire samples (Morocco and Algeria). ‘Han’ was removed because it was redundant. Hence, there were 86 populations left.

Mantel Test

The Mantel correlations for MIX-Height and SIB-Height were .291 and .395 (p < .01) respectively. The Mantel correlation for Height was .152 (p = .009).

Partial Mantel Test

A regression model (Piffer, 2015) and the partial Mantel test (Legendre & Legendre, 1998) using R package ncf (Bjornstad, 2022) were carried out to discern the effects of selection from those of ancestry.

The results of the regression model with Height (distance matrix) as dependent and PGS + F_ST (distance matrices) as independent variables are shown in Table 7.

Table 7. Regression of phenotypic distances on PGS and F_ST distances

Note: PGS, polygenic score.

***p = .001 **p = .01 *p = .05.

The partial Mantel correlation coefficient between MIX-Height and Height after controlling for F_ST was .193 (p = .003). Conversely, the partial correlation between Height and F_ST after accounting for MIX-Height was not significant and around 0 (.092).

The partial Mantel correlation coefficient between SIB-Height and Height after controlling for F_ST was .323 (p = .002). Conversely, the partial correlation between Height and F_ST after controlling for SIB-Height was around zero (.014).

Moran’s I with KNN

Global Moran’s I was calculated using the spdep R package (Bivand & Wong, 2018) employing the spatial K-nearest neighbors (KNN) method for K values ranging from 1 to 6. This KNN approach estimates the value for each case by averaging the values from the K nearest populations, as determined by the distance matrix. The results, presented in Figure 3, indicate weak to moderate spatial autocorrelation, with Moran’s I values ranging from 0.178 to 0.286.

Figure 3. Moran’s I for average body height, MIX-Height and SIB-Height for different values of K.

Spatial Autoregressive and Spatial Error Model

Spatial autoregressive and spatial error models were employed to investigate the relationship between population body height, PGSs, and spatial genetic structure. After filtering the genetic distance matrix, neighbors for each population were identified based on where the genetic distance fell below specific F_ST thresholds, resulting in a neighbour list.

For each F_ST threshold, models were constructed using average body height as the dependent variable. These models incorporated two key predictors:

1. 1. A non-spatial predictor: either MIX-Height or SIB-Height PGSs;

2. 2. A spatial component: In the spatial autoregressive model, a spatial lag term (lag Height) was computed using the F_ST weights.

In the spatial error model, the spatial error term accounted for autocorrelation in the residuals. Both models were run to capture different aspects of spatial dependence in the data. The results, presented in Table 8, showcase the p values for the PGS and the spatial autocorrelation parameters ρ (rho) and λ (lambda) for the spatial autoregressive and the spatial error models respectively. The Likelihood Ratio (LR) test evaluates whether including a spatial dependence parameter (such as ρ in a spatial autoregressive model or λ in a spatial error model) significantly improves the model’s fit compared to a model without spatial dependence. The p value associated with the LR test indicates whether the improvement in model fit (due to including the spatial parameter) is statistically significant.

Table 8. Spatial autoregressive and spatial error models with different F_ST thresholds

Note: PGS, polygenic score.

***p = .001 **p = .01 *p = .05.

Spatial Autoregressive and Spatial Error Model With HDI and PGS as Predictors

We added an environmental variable (HDI) to the spatial autoregression and error models to assess the impact of genetic and environmental factors after accounting for spatial autocorrelation. Both PGS and HDI were significant predictors in all models, whereas the autocorrelation parameter was significant only in the models with the most stringest F_ST threshold (F_ST <0.03) (Table 9).

Table 9. Spatial autoregressive and error models with HDI and PGS as predictors

Note: HDI, Human Development Index; PGS, polygenic score. Statistically significant p values are shown in bold type.

Tests of Divergent Selection

The Q_ST-F_ST test was used to identify divergent selection. We calculated Q_ST, random Q_ST (rQ_ST), and F_ST for 26 populations from the 1000 Genomes Project (Table 10). We used the Q_ST-F_ST test to detect divergent selection, calculating Q_ST, random Q_ST, and F_ST for 26 populations from the 1000 Genomes Project (Table 10). Following the notation from Le Corre and Kremer (2012), we denote the F_ST at quantitative trait loci as F_STQ and use F_ST to represent neutral genetic differentiation.

Table 10. Results of Q_ST-F_ST test

Note: GWAS, genomewide association study; MC, Monte Carlo (empirical) p value

To assess the cross-population LD component of population differentiation, we compared Q_ST to rQ_ST. The results of this simulation are shown in Figure 4.

Figure 4. Distribution of random Q_ST (rQ_ST) versus Q_ST (dashed line).

To generate rQ_ST, we shuffled effect and nFoneffect alleles 1000 times using Monte Carlo simulation. This process yielded 1000 Q_ST values. We also report the z score and p value for Q_ST/rQ_ST

As theory predicts, randomly shuffled Q_ST values closely matched the F_ST values of GWAS SNPs (.083 vs. .081 and .098 vs. .103). The Q_ST /F_ST ratio exceeded 1 for both MIX-Height and SIB-Height, but only MIX-Height showed statistical significance (p < .01).

We also carried out the local (pairwise) Q_ST test for MIX-Height. SIB-Height was not used because it failed to reach significance in the global comparison, hence it did not have the power to detect local associations after Bonferroni correction. Figure 5 shows the Q_ST values for the pairwise comparisons and the p value of the Q_ST values that were significant after bootstrapping (n = 1000).

Figure 5. Local (pairwise) Q_ST test for MIX-Height.

For the allelic differentiation component (F_ST), we compared the F_ST distribution of randomly matched SNPs to that of GWAS SNPs (Figure 6).

Figure 6. F_STQ versus F_ST (randomly matched single nucleotide polymorphisms).

Discussion

The results of this study are consistent with a model where both genetic and environmental factors influence population differences in adult body height. The results replicate across the WGS and WGS + array datasets comprising 51 and 89 (partially overlapping) populations.

This is the first study to systematically explore the impact of genetics and environmental factors on body height across multiple ethnic groups, using a comprehensive approach that integrates polygenic scores (PGS) and various environmental indicators. Previous studies have investigated the genetic determinants of body height and the role of environmental factors, but they often focused on specific populations or regions. For instance, Turchin et al. (2012) explored the genetic basis of body height in European populations, while Deaton (2007) examined the relationship between body height and socioeconomic conditions within countries. However, these studies did not combine genetic and environmental factors in a unified model across diverse global populations.

In our analysis, PGSs derived from the MIX-Height GWAS consistently demonstrated strong associations with average measured body height across both datasets. Multiple regression models using MIX-Height PGSs demonstrated significant PGS beta values, ranging from 0.402 to 0.588. These findings highlight the robust predictive power of the MIX-Height PGS in capturing genetic influences on body height.

Environmental factors such as the HDI, daily caloric intake, and protein consumption also showed significant associations with body height. For instance, HDI beta values ranged from 0.539 to 0.968, while calories and proteins had beta values ranging from 0.638 to 0.837. These results underscore the critical role of nutrition and overall living standards in determining adult body height.

The combined models incorporating both genetic and environmental factors generally yielded the highest model R² values, indicating superior explanatory power. For example, the model combining MIX-Height PGS with HDI and caloric intake achieved an R² of .846, highlighting the synergistic effects of genetics and environment on body height.

Models utilizing SIB-Height PGS also exhibited significant associations, though their PGS beta values and model R² values were generally lower compared to MIX-Height PGS models. This may reflect the smaller sample sizes and reduced phenotypic and genotypic variance within families, which limit the power of sibship GWAS to detect statistical associations.

The consistency of these results across both WGS and WGS + array datasets further strengthens our conclusions. While the WGS dataset generally produced higher model R² values, the replication of findings in the WGS + array dataset demonstrates the robustness of our models across different genomic platforms.

This study addresses the African paradox proposed by Deaton (2007), which observed high average stature in African countries despite poor nutrition and childhood health. Our findings suggest that this phenomenon is likely due to a high genetic endowment for body height among African populations, as evidenced by their higher than average polygenic scores (Figures 2a and 2b). Additionally, our results challenge Deaton’s hypotheses that (1) genetic differences do not account for substantial variation in body height among low income/high childhood mortality countries and (2) nutrition and childhood health are not significant predictors of body height variations among such countries. When we subset the dataset to only the countries with high childhood mortality, the effects of SIB-Height and MIX-Height, infant mortality and caloric intake were all significant, indicating that both environmental and genetic factors influence body height in these populations. Contrary to Deaton’s proposed balancing effect of infant mortality, we found that the effect of infant mortality on body height was negative (Table 4). When we coded HDI as a dummy variable, with a cut-off value of 0.8 separating low from high HDI populations, the interaction between MIX-Height and HDI was not significant, but the one between Calories and HDI was significantly negative, indicating, contrary to Deaton’s predictions, that the positive effect of calorie intake on body height is reduced in the higher HDI group.

Our study reveals that the Dutch population achieved the highest MIX-Height polygenic score, providing compelling genetic evidence for their status as the world’s tallest nationality. This finding lends molecular support to the theory proposed by Stulp et al. (2023) that natural selection has favored taller individuals in the Netherlands over time.

Interestingly, our analysis of nutritional factors yielded unexpected results. Contrary to popular belief, the Dutch protein and calorie intake levels were found to be average when compared to other developed countries. This observation challenges the widespread notion that the exceptional body height of the Dutch population is primarily attributable to high consumption of dairy products.

Our models incorporated latitude as a climate proxy, drawing on Bergmann’s rule, which posits that colder environments lead to larger body sizes, including increased body height in humans. While the global correlation was weak, our analysis revealed a more pronounced effect of latitude on body height PGSs when accounting for major ancestry groups (‘superpopulations’). This effect reached statistical significance in three out of four models (Table 6). Notably, African ancestry showed a strong positive influence on body height PGS, suggesting that African populations are taller than expected given their relatively equatorial geographical origin. This finding positions African populations as an outlier in the latitude-body height relationship. Conversely, East Asian populations appear to have a genetic predisposition for shorter stature than their latitude would predict, though this effect is generally less pronounced than the African outlier. These unexpected results invite further exploration of potential underlying mechanisms, such as sexual selection pressures or adaptations for disease resistance. It is worth noting that Allen’s rule (Allen, 1877), which proposes that taller, leaner body types can be advantageous for heat dissipation in hot climates, offers predictions that contrast with Bergmann’s rule. Specifically, Allen’s rule would predict that people living in hotter climates have relatively longer limb length, but not necessarily taller stature.

This complexity underscores the need for a more nuanced understanding of the interplay between genetics, climate and human body height variation.

Our study incorporated spatial autocorrelation — more accurately termed ‘genetic autocorrelation’ in this context — using several methodologies: Mantel tests, partial Mantel correlations, Moran’s I statistics, spatial autoregressive and spatial error models. The Mantel tests (Figures 1a and 1b) and Moran’s I analyses (Figure 3) revealed evidence of moderate genetic autocorrelation. However, when we integrated genetic autocorrelation into regression models alongside PGS as predictors, its effect was diminished. Specifically, in Mantel regression and partial correlation analyses, the impact of genetic autocorrelation was minimal (Table 7).

Interestingly, phenotypic height exhibited much stronger spatial autocorrelation than genotypic height (i.e., PGS), as shown by the Moran’s I values for different values of k in Figure 3. This was despite ‘spatial’ autocorrelation being actually measured from genetic distance matrices, thus not directly indicating geographical or cultural proximity.

Possible explanations for this phenomenon is that countries that are genetically similar (as indicated by F_ST distances) often share similar environments due to geographical proximity, historical migration patterns, or socioeconomic similarities. These shared environmental factors could amplify genetic predispositions, leading to more similar phenotypic heights across populations, even if PGS values for height are not perfectly aligned. There could be a correlation between genetics and the environment within genetically similar populations. For example, countries with a shared genetic background (as measured by F_ST) may also have similar cultural or socioeconomic conditions that influence height. This gene-environment correlation can cause genetically similar countries to exhibit more phenotypic similarity than would be predicted by their PGS alone. In this case, the environment reinforces the genetic potential for height, leading to stronger autocorrelation in phenotypic height.

In our regression analysis, we opted to define neighbors based on F_ST distance rather than using a fixed number of neighbors, as done in the KNN method, which can include ‘neighbors’ with large genetic distances. This choice was supported by the Mantel correlograms (Figures 1a and 1b), which revealed negative autocorrelation at higher F_ST values (>.05), potentially reducing our ability to detect autocorrelation.

The spatial autoregressive and error models showed a stronger autocorrelation effect when using MIX-Height compared to SIB-Height as the predictor, reflected in rho and lambda coefficients of approximately 0.7 (Table 8). Despite this, PGS remained a significant predictor. The standardized beta for MIX-Height ranged from 0.279 (p = .004) in the spatial autoregressive model with F_ST < 0.03, to 0.435 (p < .001) in the spatial error model with F_ST < 0.09.

For SIB-Height, PGS was a significant predictor across all F_ST thresholds, with beta values ranging from 0.328 in the spatial autoregressive and error models with F_ST < 0.03, to 0.783 in the spatial error model with F_ST < 0.07. However, the genetic autocorrelation effect was significant only at the lowest F_ST threshold of .03.

The degree of autocorrelation decreased after incorporating an environmental predictor (HDI) into the regression models, indicating that phenotypic-level autocorrelation is partially driven by environmental factors. In models with the lowest F_ST threshold, rho decreased from 0.698 to 0.422 for MIX-Height and from 0.588 to 0.373 for SIB-Height (Table 9). Notably, PGS and HDI remained significant predictors of average height in all models that accounted for spatial autocorrelation across all F_ST thresholds

From a theoretical evolutionary genetics perspective, controlling for genetic autocorrelation serves as a method to detect the overdispersion of genetic values among populations relative to neutral expectations, as reflected in the F_ST distances matrix. Significant levels of overdispersion are interpreted as signals of local adaptation. Conceptually, this approach bears a strong resemblance to Q_ST-F_ST comparisons, a widely employed test for divergent selection (Merilä & Crnokrak, 2001).

Q_ST values exceeded Fst(q) for both MIX-Height and SIB-Height, with Q_ST-F_ST ratios of 3.41 and 1.67 respectively. A Monte Carlo simulation, involving reshuffling of effect and noneffect alleles, showed significant results for MIX-Height (p = .008) but not for SIB-Height (p = .179) (Figure 5).

This test alone is overly conservative because selection effects encompass both allelic differentiation and the covariance of allelic effects across populations. Allelic differentiation influences the difference between F_ST at GWAS loci (F_STQ) and F_ST at neutral loci, while Q_ST randomization reveals the impact of allelic effect covariance across populations — the primary component of differentiation in polygenic traits (Berg & Coop, 2014; Le Corre & Kremer, 2012;).

To address this, we conducted the F_ST enrichment (F_STQ-F_ST) test to examine the allelic differentiation component. The results were significant for MIX-Height but not for SIB-Height (p < .001 and .10, respectively) (Figure 6).

It is crucial to note that the interpretation of PGSs within a causal framework is only valid within this specific theoretical context. This caveat is particularly salient for PGSs that account for only a small proportion of phenotypic variation, such as those derived from sibship GWASs. The underlying reason for this limitation lies in the complex interplay between local adaptation and allele frequencies across populations.

When local adaptation occurs, it tends to increase the frequencies of alleles that have similar effects (either increasing or decreasing) on a given trait across populations. Importantly, this phenomenon is not limited to the alleles identified by GWAS; it extends to other unidentified alleles that contribute to the trait. Consequently, the effect of the genetic variants comprising the polygenic score on the population-level trait is not direct. Instead, it is mediated by a broader set of alleles that, under evolutionary expectations, follow similar patterns of frequency distribution (Berg & Coop, 2014; Piffer, 2013).

The underlying reason for this limitation lies in the complex interplay between local adaptation and allele frequencies across populations.

This indirect relationship can be further elucidated as follows:

GWAS identifies a subset of alleles associated with a trait, which are used to construct the PGS. These identified alleles represent only a fraction of the total genetic variation influencing the trait.

Local adaptation influences the frequencies of these identified alleles across populations. This process occurs through natural selection acting on the phenotypic effects of these alleles in different environmental contexts.

Simultaneously, local adaptation affects the frequencies of unidentified alleles that also contribute to the trait. These unidentified alleles, while not captured in the PGS, are subject to the same evolutionary pressures as the identified alleles.

The observed effect of the PGS at the population level is thus a proxy for the cumulative effect of both identified and unidentified alleles. This relationship arises because the PGS, based on a limited set of identified alleles, serves as an indicator of broader genetic patterns shaped by local adaptation.

This cumulative effect reflects the broader evolutionary forces shaping genetic variation related to the trait, rather than a direct causal relationship between the specific alleles in the PGS and the trait variation across populations. In essence, the PGS acts as a marker for the overall genetic architecture influenced by local adaptation, encompassing both the measured and unmeasured genetic components.

Understanding this indirect relationship is crucial for accurately interpreting the results of studies employing PGSs, particularly in the context of population genetics and local adaptation. It underscores the need for caution when making causal inferences based on PGSs, especially when these scores explain only a small fraction of phenotypic variation. Researchers must consider the broader evolutionary context and the potential influence of unidentified genetic factors when drawing conclusions from PGS analyzes across populations.

In summary, our study provides compelling evidence that both genetic predispositions and environmental conditions significantly contribute to adult body height. The integration of PGSs and environmental factors offers a comprehensive understanding of the determinants of body height, reinforcing the importance of considering both genetic and nongenetic factors in future research on human growth and development. However, our study had several limitations:

Population diversity. Although the study includes a wide range of populations, the sample sizes for certain regions may be insufficient to capture the full genetic and environmental variability within those areas. This could lead to biased estimates of the impact of both genetic and environmental factors on body height.

Environmental measures. The environmental indicators used in this study, such as HDI, daily caloric intake and protein consumption, are proxies and may not fully capture the complex environmental influences on body height. More granular data on specific health and nutritional factors would provide a clearer picture.

Model limitations. The models used in this study, while comprehensive, may still be overly simplistic in capturing the interplay between genetics and environment. Nonlinear relationships and interactions between different environmental factors could be more thoroughly explored in future research.

Genetic data representation. The genetic data, although extensive, may not cover all relevant genetic variants influencing body height. Future studies could benefit from even more comprehensive genomic data, including rare variants and epigenetic factors. Moreover, the use of F_ST distance matrices, albeit common in population genetics, is not entirely correct as it does not always satisfy the triangle inequality and thus is not a metric (Arbisser & Rosenberg, 2020).

Replication and validation. While the results are robust across the datasets used, replication in other independent datasets and validation in different demographic contexts would strengthen the generalizability of the findings.

Ethical considerations. The study of genetic differences across populations must be conducted with careful consideration of ethical implications and potential misuse of the findings. This includes ensuring that the research does not reinforce stereotypes or lead to discrimination.

Accuracy of measured height. The accuracy of average measured body height data can vary, and discrepancies in measurement techniques across different studies or populations might affect the results. Consistency in measurement methods is crucial for reliable comparisons.

Missing data. Several ethnic groups had missing data on average body height, which could skew the results. The absence of complete data for all populations limits the study’s ability to fully capture the genetic and environmental determinants of body height across diverse groups.

These limitations should be considered when interpreting the results and in designing future research to build on these findings.