Adaptive diversification of growth allometry in the plant Arabidopsis thaliana.

Seed plants vary tremendously in size and morphology; however, variation and covariation in plant traits may be governed, at least in part, by universal biophysical laws and biological constants. Metabolic scaling theory (MST) posits that whole-organismal metabolism and growth rate are under stabilizing selection that minimizes the scaling of hydrodynamic resistance and maximizes the scaling of resource uptake. This constrains variation in physiological traits and in the rate of biomass accumulation, so that they can be expressed as mathematical functions of plant size with near-constant allometric scaling exponents across species. However, the observed variation in scaling exponents calls into question the evolutionary drivers and the universality of allometric equations. We have measured growth scaling and fitness traits of 451 Arabidopsis thaliana accessions with sequenced genomes. Variation among accessions around the scaling exponent predicted by MST was correlated with relative growth rate, seed production, and stress resistance. Genomic analyses indicate that growth allometry is affected by many genes associated with local climate and abiotic stress response. The gene with the strongest effect, PUB4, has molecular signatures of balancing selection, suggesting that intraspecific variation in growth scaling is maintained by opposing selection on the trade-off between seed production and abiotic stress resistance. Our findings suggest that variation in allometry contributes to local adaptation to contrasting environments. Our results help reconcile past debates on the origin of allometric scaling in biology and begin to link adaptive variation in allometric scaling to specific genes.

At the core of the quest for understanding and predicting biological diversity is the apparent paradox that, despite the phenotypic changes that underlie divergent ecological strategies, there seem to be constant or near-constant parameters across life forms (1). The latter is assumed to result in part from biophysical constraints limiting the range of possible trait values (2), as well as from strong stabilizing selection for optimal phenotypes (3, 4). Consistently, body size variation in multicellular organisms is associated with many scaling regularities. Kleiber (5) first reported that the consumption of energy (metabolic rate, G) varies to the 3/4-power of organism mass M, such that G = GM, implying that a 10-fold increase in M produces a 5.6-fold increase in G in virtually all organisms.

Several physiological models have been proposed to explain this constancy. The most prominent of these is metabolic scaling theory (MST) (6), which predicts that scaling exponents of several traits tend to take on “quarter-power” values (e.g., 3/4, 1/4) as the outcome of an optimal balance between the scaling of hydraulic transport costs and the scaling of exchange surface areas (e.g., leaf area in plants) (7). According to MST, the scaling of physiological rates matches the ability of exchange surfaces to obtain resources from the environment and then distribute them to metabolizing cells through the vascular network. Because the branching geometry of this network is highly constrained in space, it is predicted that selection that minimizes the costs of resource transport and at the same time maximizes the uptake of resources will lead to “allometrically ideal” organisms characterized by a common set of quarter-power scaling relationships with body mass.

Empirical observations support MST predictions across land plants, where several traits, including organismal growth rate, scale as body mass raised to the power of 3/4 (8, 9). On the other hand, the scaling exponent can vary across plants (10–12), or scaling can be constant but can deviate from 3/4 (13). These seemingly contradictory observations have been proposed to reflect (i) phenotypic differences, like life history and differences between species or populations (9, 10); (ii) physiological changes along environmental gradients (14, 15); or (iii) nonlinearity in hydrodynamic resistance and metabolic scaling (16). Thus, important questions about the evolution of allometry remain (4). For example, is the prevalence of ubiquitous scaling relationships the result of stabilizing selection acting to remove unfit genetic allometric variants? Does variation in the scaling exponent reflect adaptation and genetic diversification, or developmental plasticity?

To address these and related questions, we examined how growth rate scales with body size in a genetically diverse population of Arabidopsis thaliana accessions (Dataset S1), a species that exhibits three orders of magnitude in plant dry mass (10) and occurs in a wide range of contrasting environments (17). We provide evidence indicating that scaling variation is maintained by an adaptive trade-off between alternative environments. We show that this variation has a polygenic basis, and that there is genetic correlation between allometry and local climate.

Results

Variations of A. thaliana Growth Scaling with Climate.

The scaling exponent of growth is conventionally quantified as the slope, θ, of the allometric function y = α + θx, where x and y are the logarithms of plant biomass and absolute growth rate (GR), respectively. Fitting the allometry of the mean absolute GR (in milligrams per day), estimated as the ratio of final plant dry mass (in milligrams) over duration of the life cycle (in days), across A. thaliana accessions returned a scaling exponent, θ, that is not significantly different from the MST predicted value of 3/4 (y = −1.07 + 0.74x; r2 = 0.97; slope CI95% = [0.725, 0.750]) (Fig. 1). This value is the same as that observed across vascular plant species (box in Fig. 1); however, the relationship is not a pure power function, and instead is better explained by a nonlinear quadratic function [y = −1.93 + 1.43x − 0.14x2; ΔAIC (Akaike information criterion) = −192.4], (Fig. 1). Our analyses indicate that this curvilinear scaling relationship was due to differences in θ between accessions, which can be estimated as the first derivative of the quadratic function (θ = 1.43 − 0.27x), and which varied among accessions from 0.47 to 1.10 (Fig. 1 and Fig. S1). The broad-sense heritability, H2, of θ was 0.95, which is higher than any other trait measured in this study (Table S1), indicative of a high amount of variance explained by genetic effects in our highly controlled growth conditions.

Fig. 1.

Variation of growth scaling in A. thaliana. (A) Linear (dashed line) and quadratic (solid line) fits of mean GR vs. final dry mass in 451 A. thaliana accessions. (Box) Linear fit (black line) of GR vs. plant dry mass in 333 vascular plant species from Niklas and Enquist (8). (B) Distribution of the scaling exponent derived from the quadratic fit in the 451 A. thaliana accessions. (C) Relationship among RGR at growth maximum, plant lifespan, and scaling exponent in the 451 accessions. The black curve represents loess fit ±95% CI (gray area). In all panels, dots and triangles represent genotypic and species means, respectively, color-coded by the value of the scaling exponent reported in B.

Modeling the dynamics of plant dry mass accumulation from imaging data (18) revealed that the estimated relative growth rate (RGR) explains 18% of the variation in the scaling exponent (P < 0.001), with both negatively correlated with plant lifespan (P < 0.001) (Fig. 1 and Dataset S2). Previous studies have shown that variation in A. thaliana growth allometry is positively correlated with the carbon assimilation rate and nutrient concentration, but negatively correlated with lifespan (10). Thus, variation of growth allometry in A. thaliana connects life history variation to the strategies for leaf resource use. At one end of the distribution are high scaling exponents, representative of “live fast/die young” strategies that maximize resource capture (i.e., high RGR and carbon assimilation rate) at the expense of plant lifespan and final size. At the other end are low scaling exponents, representative of “live slow/die old” strategies that maximize the retention (i.e., thick leaves with low nutrient concentration and long lifespan) rather than the acquisition of resources.

We then examined the correlations between the scaling exponent and 21 climatic variables, which include 19 Bioclim variables (www.worldclim.org/bioclim), as well as the estimated mean annual potential evapotranspiration (PET, in millimeters) and the aridity index (19) at the geographic origin of the accessions. Consistent with the idea that resource-acquisitive plants (i.e., early-flowering/fast-growing ecotypes) are more adapted to hotter and drier regions, the scaling exponent was positively correlated with the mean annual temperature measured at the collection point of the accessions (Fig. 2 and Dataset S2). The strongest correlations were with the maximum temperature of the warmest month and the mean temperature of the warmest quarter (r = 0.30 and 0.28, respectively) (Fig. 2 and Dataset S2). Inversely, the scaling exponent was negatively correlated with precipitation, specifically with precipitation during the driest quarter (Fig. 2), precipitation seasonality, and the aridity index (Dataset S2). In contrast, it was not correlated with the altitude at the collection point.

Fig. 2.

Relationship between scaling exponent and climate. (A–C) Correlations between the scaling exponent measured across the 451 accessions and local mean annual temperature (A), maximum temperature of the warmest month (B), and precipitation of the driest month (C). Dots represent the genotypic mean. Fitted lines are SMA regressions. r is the Pearson coefficient of correlation with associated P value. (D) Geographic distribution of the scaling exponent across Europe in A. thaliana, modeled as a function of 13 Bioclim variables. Colors indicate the predicted value of the scaling exponent. Black dots represent geographic origins of the accessions phenotyped.

Using stepwise regression, we found that 13 climatic variables explained >27% of the allometric variation. Four of these variables are related to summer climate and two are related to winter climate. The strongest effects were estimated for annual mean temperature, isothermality, and mean summer temperature. Modeling the geographic distribution of scaling exponent with the 13 top-correlated climatic variables as predictors showed that intermediate exponents are more common in temperate regions (Fig. 2), while extreme exponents are favored under more stressful conditions (e.g., high altitude, high latitude).

Fitness Costs and Benefits of Allometric Variation.

The scaling exponent was correlated with resource use traits, including RGR and lifespan, as well as performance-related traits, such as fruit number, a proxy for lifetime fitness in annual species. Fruit number varied from 18 to 336 per plant (Table S1 and ). However, the relationship between fitness and the scaling exponent under nonlimiting Raspberry Pi Automated Phenotyping Array (RAPA) conditions was not linear (Fig. 3). Instead, fruit number was a bell-shaped function of the scaling exponent, peaking for plants with an exponent around 3/4 and declining toward higher or lower exponents. Thus, in A. thaliana, genetic deviations from the 3/4 scaling exponent are associated with extreme resource use strategies and a general decline in fruit number (r = −0.62; P < 0.001) (Dataset S2). A polynomial regression of relative fitness using fruit number (standardized by the population mean over the scaling exponent) returned a significant negative second-order coefficient (y = 1.00 + 4.23x − 4.06x2; P < 0.001 for all coefficients), i.e., an estimate of quadratic selection gradient, |γ|, that might be indicative of stabilizing selection for the allometric exponent under benign conditions (20).

Fig. 3.

Relationships among the scaling exponent, fitness, and resistance to abiotic stress. (A) Relationship between fruit production and scaling exponent in the 451 accessions. The black curve represents loess fit ±95% CI (dashed lines). (B) Stress resistance expressed as the log10 value of the ratio of final rosette dry mass under water deficit (WD), high temperature (HT), and both (WDxHT) compared with control (cont) conditions, across 120 A. thaliana recombinant inbred lines. The data have been published previously (10, 23). Dots indicate genotypic means (n = 4). Colored curves represent loess fit ±95% CI (dashed lines).

Conversely, deviation from 3/4 scaling was positively correlated with survival under severe drought (r = 0.16, P < 0.05; measured in ref. 21 across 210 common accessions) (Dataset S2) and negatively correlated with growth reduction under moderate drought (r = −0.26, P < 0.05; measured in ref. 22 across 60 common accessions) (Dataset S2). However, neither stress-resistance trait was correlated with the scaling exponent itself. This suggests that deviation of allometric exponents from 3/4 in either direction is associated with increased resistance to stressful conditions at the expense of reduced reproductive fitness under benign conditions.

Consistently, reanalysis of an experimental population phenotyped for tolerance to combined high temperature and water deficit (23) pointed to higher stress sensitivity of accessions with scaling exponents close to 3/4 (Fig. 3). In contrast, allometric exponents at both the low end and high end of the distribution were correlated with improved stress tolerance, specifically under high temperature (Fig. 3). A possible explanation of this result could be that a “fast” strategy with high scaling exponents allows stress escape by maximizing resource acquisition and completion of the life cycle before a short window of nonstressful conditions closes (23). Alternatively, the “slow” strategy might support stress tolerance by reducing metabolic activities and thus the resource demand associated with fast growth (10).

Genetic and Evolutionary Bases of Allometric Variation.

Because we suspected that allometric variation might result from adaptation to the diverse environments at the places of origin of accessions, we looked for genetic evidence of local adaptation and of genetic diversification with climate. Principal component analysis (PCA) performed after eigendecomposition of the relatedness matrix revealed that the scaling exponent was correlated with population structure, notably with the second PCA axis (r = 0.37; P < 0.001), which explains 28% of the total genetic variation and mainly differentiates accessions from the ancestral (“relict”), N. Sweden, and Spain genetic groups (17) (Fig. S2). In contrast, flowering time was correlated with the first PCA axis, which explains 42% of genetic variation and is associated with longitudinal divergence among accessions (Fig. S2). Compared with the relict genetic group (17), the scaling exponent differed significantly (P < 0.001) for the N. Sweden and S. Sweden groups, but not for the other eight groups (P > 0.3). The Qst of the scaling exponent, measured as the ratio of between-group phenotypic variance to total variance, was above the 0.9 quantile of genome-wide Fst (Qst/Fst ratio = 2.14; P < 0.001) (Fig. S3 and Table S1), which is potentially indicative of polygenic selection acting on the scaling of plant growth (24).

We ran genome-wide association (GWA) models on the scaling exponent θ and the 21 climatic variables using EMMAX to correct for population structure (25). After multiple-testing correction, a total of 8,250 single nucleotide polymorphisms (SNPs) out of 1,793,606 tested were significantly associated with at least one phenotypic trait or climatic variable (26) (Dataset S3). Only six SNPs were significantly associated with the scaling exponent [false discovery rate (FDR) <0.05]. Five of these six SNPs were located in the same region on chromosome 2 (Fig. 4) and were associated with the maximum temperature of the warmest month (Fig. 4). Three SNPs were also significantly associated with the mean annual temperature and the mean temperature of the coldest month (Dataset S4). The same genomic region showed a strong association with precipitation during the driest month (Fig. 4), although the six SNPs associated with scaling variation did not reach the significance threshold for this climatic variable (FDR >0.05). In contrast, no significant SNPs were shared between RGR, lifespan, fruit number or rosette dry mass, and the climatic variables (Dataset S4), suggesting that a genetic association between traits and climate is relatively rare.

Fig. 4.

GWA mapping of allometric variation in A. thaliana. (A–C) Test statistics for SNP associations (EMMAX) with the scaling exponent (A), maximum temperature during the warmest month (B), and precipitation during the driest month (C). Dots represent the top 1% associated SNPs along the five chromosomes (alternate gray and black dots represent chromosomes). Orange lines represent the genome-wide significance threshold with Bonferroni correction at α = 0.05 (solid line) and α = 0.1 (dashed line). The red triangle is PUB4 (FDR <0.05). (D and E) Correlation between SNP effects (BSLMM) for the scaling exponent and maximum temperature of the warmest month (D) and precipitation of the driest month (E). Black dots represent a similar SNP effect for x and y variables (both positive or both negative). r is Pearson’s coefficient of correlation (***P < 0.001).

One SNP among the five associated with both the scaling exponent and the maximum temperature of the warmest month was located in the U-box protein gene PUB4 [At2g23140; minor allele frequency (MAF) = 6.1%] (Fig. 4). As E3 ubiquitin ligases, U-box proteins are involved in protein turnover, a key regulatory component of plant responses to abiotic stresses (27). PUB4 plays a notable role in a quality control pathway that removes damaged chloroplasts (28). Two other SNPs were located in the nearby cytochrome P450 gene CYP81D6 (At2g23220), 40 kb from PUB4 (r2 = 0.63). CYP450s catalyze the production of diverse secondary metabolites involved in biotic and abiotic stress responses (29). The remaining two SNPs were also linked to PUB4 and CYP81D6 but affected noncoding sequences.

We note that the PUB4 polymorphisms account for only approximately 1% of the genetic variance in the scaling exponent. Because broad-sense heritability was H2 >95%, many other loci are expected to contribute to allometric variation, potentially reducing the power of classical GWA to detect SNPs significantly associated with the scaling exponent. For instance, given the strong correlation between the scaling exponent and plant lifespan (Dataset S2), we expected to find that many flowering time genes would be significantly associated with allometry; however, no SNP reached the significance threshold for lifespan in our analysis (FDR >0.05), and thus we do not have evidence identifying flowering time genes as predictive of allometric variation. This might be due to overcorrection for population structure or to the high number of SNPs involved in phenotypic variation among accessions. Indeed, a strong correction for population structure might be inappropriate if many genes across the entire genome contribute to the phenotype in question.

To account for the potentially complex genetic architecture of traits, we ran Bayesian sparse linear mixed models (BSLMMs) implemented in GEMMA (30). BSLMM models two hyperparameters: a basal effect, αi, that captures the fact that many SNPs contribute to the phenotype, and an extra effect, βi, that captures the fact that not all SNPs contribute equally. SNP effects, which can be estimated as the sum of αi and βi (30), were strongly correlated between the scaling exponent and all climatic variables except annual temperature range (Dataset S5). As expected, correlations between SNP effects on scaling exponent and climate were strongest for mean annual temperature and temperature and precipitation during summer (Dataset S5). Consistent with the measurement of broad-sense heritability (H2), “chip heritability” (a proxy for narrow-sense heritability, h2, measured with GWA) was very high for the scaling exponent (h2 = 0.87 vs. H2 = 0.95) (Table S1), suggesting that most of the phenotypic variance can be explained by the additive effects of SNPs controlling allometric variation.

Gene Ontology (GO) analysis (31) of the top 1% of genes affecting the scaling exponent revealed enrichment in genes with catalytic activity and genes related to carbohydrate metabolism, postembryonic development, posttranslational protein modification, and response to abiotic stimulus (Fig. S4 ). A large fraction of the proteins encoded by these genes are predicted to localize to plasma membranes or the chloroplast (Fig. S4). Fst values across the top-100 genes were significantly higher than genome-wide Fst values (0.23 vs. 0.17; P < 0.001) (Fig. S3), which is consistent with the Qst analysis and indicative of polygenic selection on the genes controlling growth allometry. As expected, PUB4 is among the top-100 genes associated with plant allometry, showing strong effects on both the scaling exponent and climatic variables (Fig. 4 ). We estimated that PUB4 alone favors plant adaptation to warmer and drier summers by up to +1.4 °C and −3 mm (Fig. 4 ) through an increase of the scaling exponent by up to +0.03.

A scan for genomic signatures of selection in the 50-kb region around PUB4 revealed increased Tajima’s D (Fig. 5) and SNP-level Fst (Fig. 5) values, but no signatures of recent selection sweeps. Tajima’s D is an index of allelic diversity that quantifies departures from the standard neutral model (32). High values indicate an excess of intermediate-frequency alleles, a potential sign for balancing selection, specifically in A. thaliana, in which Tajima’s D is commonly negative due to recent population expansion and selfing (33, 34). This finding is consistent with molecular signatures of climate selection previously observed in A. thaliana (35, 36). Moreover, climate envelope modeling of PUB4 allelic distribution revealed a strong geographic structure associated with summer conditions; the major PUB4 allele is found mostly in temperate and cold northern parts of Europe (Fig. 5), while the minor allele is mostly Mediterranean (Fig. 5). This supports the role of PUB4 in evolutionary adaptation to warmer and drier regions around the Mediterranean through variation in growth scaling.

Fig. 5.

Genomic signatures of adaptation to climatic conditions at genes controlling the scaling exponent. (A and B) Tajima’s D (A) and Fst (B) in a 50-kb region around PUB4 and CYP81D6. Gray dots are mean values in 1-kb bins; red lines indicate positions of significant SNPs. (C and D) Predicted geographic frequency of the major (C) and minor (D) alleles at PUB4 following climate-envelope modeling with 19 Bioclim variables. The color gradient indicates predicted allele frequency.

Discussion

Metabolic allometry links physiology, ecology, and evolution at different levels of organization (4, 6, 37, 38). The study of scaling relationships in both plants and animals is grounded on the importance of universal metabolic properties that allow the measurement and prediction of critical rates of energy flow from individuals to the biosphere (6, 39). Nonetheless, explanations of the origin of allometric variation between species remain elusive, despite the recognized role of evolutionary processes in animals (40). Changes in the allometric scaling intercept G0 in response to selection are well documented (41), but evidence of the evolution of allometric slopes θ is scarce (but see ref. 42), particularly in plants, where the focus has been on the specific value that the allometric slope should take (e.g., 2/3 vs. 3/4 vs. 1) (9, 13, 43).

Our present findings reconcile the recent debate on the origin of biological allometry. On the one hand, our results support the ideas that growth allometry varies significantly, and that genetic variation in allometry is maintained within species. On the other hand, the canonical 3/4 scaling exponent reported within and across plant and animal species was found to be associated with a phenotypic optimum that maximizes fitness under benign conditions, consistent with a role of stabilizing selection (4). Nonetheless, depending on the local environment, deviations in both directions from the 3/4 scaling exponent might be advantageous for stress resistance despite their cost to seed production. Thus, stabilizing selection on metabolic allometry could be disruptive under unfavorable environments, as we have found for A. thaliana. Allometric adaptation may be due to, for instance, selection for fast growth and a short lifespan to escape drought, or selection for resistance to hydraulic cavitation associated with reduced stomatal conductance and carbon assimilation in late-flowering ecotypes (23, 44).

Specifically, these findings provide insight into the important role of allometry for local adaptation to various climates in A. thaliana. Moreover, our results inform our understanding of the evolutionary basis of the tenets of MST. The maintenance of high intermediate-frequency nucleotide diversity in genes affecting allometry could result from long-term, geographically heterogeneous selection to optimize growth and survival in contrasting environments. This appears to have resulted in the genetic diversification of the scaling exponent around the intraspecific and interspecific mean of 3/4, potentially reconciling the original MST prediction of an optimal scaling 3/4 value with observed departures from it that have generated past debates (45).

An intriguing question is whether the observed variation in scaling exponents across species (46) is associated with a similar climate adaptation as we observed for A. thaliana. Interspecific and intraspecific variation in the vascular network and its impact on hydrodynamic resistance, resource distribution, and plant allometry is already being explored (47, 48). If genetic variability in growth allometry is confirmed in other species and associated with the environment (e.g., ref. 11), this would have important implications for our understanding of the physiological bases of plant adaptation. Moreover, it would connect macroevolutionary patterns of trait covariation observed across species to microevolutionary processes occurring within species.

Materials and Methods

Published Data.

For stress resistance analysis, we used published data from two studies on the response of A. thaliana natural accessions to drought, including one study in which 210 accessions shared with our study were subjected to severe, lethal drought and survival was estimated for all accessions (21), and another in which 60 shared accessions were subjected to 7 d of nonlethal drought and fresh weight was measured (22). We also reanalyzed previously published phenotypic data (10, 23) from a population of 120 Ler-2 × Cvi recombinant inbred lines (RILs) (49) grown under water deficit and high temperature.

Climatic data consisted of 19 bioclimatic variables (www.worldclim.org/bioclim) with a 2.5 arc-min resolution for the 1950–2000 CE period, plus the mean annual PET (in millimeters) and annual aridity index downloaded from www.cgiar-csi.org/data/global-aridity-and-pet-database (19). Monthly averages were calculated with 30 arc-s (ca. 1 km). Additional details are provided in .

Plant Material and Growth.

We selected 451 natural accessions of A. thaliana from the 1001 Genomes Project (17) (1001genomes.org/) (Dataset S1). Seeds were obtained from parents propagated under similar conditions in the greenhouse. Four replicates of each accession were grown, with one replicate each sown on four consecutive days. Two replicates per accession were harvested as 16-d-old seedlings for dissection, imaging, and weighing, and two other replicates were cultivated until the end of the life cycle (i.e., until fruit ripening) for trait measurement. Plants were cultivated in hydroponic culture on rockwool. Seedlings were vernalized for 4 °C (8-h light) for 41 d and then transferred to 16 °C (12-h light). Additional details are provided in .

Plant Measurements.

The RAPA system was used for continuous imaging using 192 microcameras (OV5647; OmniVision), which simultaneously acquired six daily top-view 5-megapixel images for each tray of 30 plants during the first 25 d after vernalization. Recording and storage of images were managed through embedded computers (Raspberry Pi rev. 1.2; Raspberry Pi Foundation). Inflorescences and rosettes of mature plants were separated and photographed with a Canon EOS-1 camera. The rosette was dried for at least 3 d at 65 °C and then weighed with a microbalance (XA52/2X; Rauch).

Plant lifespan was measured as the duration between the appearance of the two first leaves and the end of reproduction (additional details in ). Fruits (siliques) were counted by eye on inflorescence images of 352 plants harvested at maturity. We analyzed the inflorescence pictures of all harvested plants with ImageJ (50) to estimate the number of fruits through image 2D skeletonization (18). The inferred variables were used to predict fruit number with linear regression (glm) tested on the 352 plants for which we had both measurements (18).

Drought survival index values were obtained from published data, measured as the quadratic coefficient of the polynomial regression between green leaves and time after the end of watering; more negative values indicate lower survival (21). Measurements of growth reduction under moderate drought were also obtained from published data, measured as the percentage of rosette fresh weight after 7 d of water deficit compared with control (22). In the reanalysis of the population of 120 RILs previously phenotyped for the growth scaling exponent (10) and trait plasticity in response to a water deficit and high temperature (23), we measured the resistance to combined stresses through the log ratio of dry mass under stress or no stress. Additional details are provided in .

Modeling Growth and RGR.

Absolute GR (in milligrams per day) was estimated as the ratio of final rosette dry mass to plant lifespan. Using rosette dry mass estimated from image analysis (18), we fitted a sigmoid curve as a three-parameter logistic equation (51) with the function nls in R. From the parameters of the fitted function of each individual, we measured RGR (rosette growth rate divided by rosette dry mass, in milligrams per day per gram) at the inflection point of the growth trajectory (18).

Statistical Analyses.

All statistical analyses except genomic analyses were performed in R (52). The coefficients of correlation (and their associated P values) reported between phenotypic traits and climatic variables were the Pearson’s product moment coefficients obtained with the function cor.test in R. The effect of population structure on the scaling exponent was tested with ANOVA, using the nine genetic groups identified in the 1001 Genomes Project dataset (1001genomes.github.io/admixture-map/) after removing admixed accessions (17). Broad-sense heritability (H2) was measured as the proportion of variance explained by genotype (Vg) over total variance (Vg + Ve) in a linear mixed model fitted with the lme4 R package, such as H2 = Vg/(Vg + Ve). Similarly, Qst was measured as the amount of variance in phenotypes explained by genetic group membership. As for H2, we used linear mixed model in the lme4 package in R to fit traits against genetic groups (nine genetic groups after removal of admixed accessions).

Genetic Analyses.

Conventional GWA studies were performed with easyGWAS (25) (https://easygwas.ethz.ch/). We used 1,793,606 SNPs with an MAF >0.05 to compute the realized relationship kernel from the full sequence of the accessions (1001genomes.org/). Association analyses were performed with EMMAX (53). For polygenic GWA, we used the BSLMM implemented in GEMMA (30). GO analysis was performed online using AgriGO (bioinfo.cau.edu.cn/agriGO/) (31) and REVIGO (revigo.irb.hr/) (54).

Before Fst calculation, genetic groups in the 1001 Genomes Project collection were defined by ADMIXTURE clustering (17, 55) (1001genomes.github.io/admixture-map/). Genome-wide estimates of Weir and Cockerham’s Fst (56) were obtained with PLINK v1.9 (57). Local selection scans (Tajima’s D and Fst) were obtained in 1-kb sliding windows in the 50-kb region around PUB4 using PLINK. Selection sweep scans were carried out using SweeD software (58). Additional details are provided in .

Modeling Geographic Distribution.

We performed stepwise regression to identify the set of climatic variables that best explain the variation of the scaling exponent between 36°N and 64°N and between 10.5°W and 27.5°E. We then used linear regression of the scaling exponent with the 13 best climatic variables to predict the exponent at each location across Europe. Geographic representation was obtained with the raster package in R. We performed climate-envelope modeling of allelic frequency at PUB4 with maxent modeling (59), using the dismo and raster packages in R. We used the 19 Bioclim variables downloaded from Worldclim database at the origin of accessions with a 2.5-arc-min resolution. Additional details are provided in .