Functional trait divergence and trait plasticity confer polyploid advantage in heterogeneous environments

Summary Polyploidy, or whole‐genome duplication often with hybridization, is common in eukaryotes and is thought to drive ecological and evolutionary success, especially in plants. The mechanisms of polyploid success in ecologically relevant contexts, however, remain largely unknown. We conducted an extensive test of functional trait divergence and plasticity in conferring polyploid fitness advantage in heterogeneous environments, by growing clonal replicates of a worldwide genotype collection of six allopolyploid and five diploid wild strawberry (Fragaria) taxa in three climatically different common gardens. Among leaf functional traits, we detected divergence in trait means but not plasticities between polyploids and diploids, suggesting that increased genomic redundancy in polyploids does not necessarily translate into greater trait plasticity in response to environmental change. Across the heterogeneous garden environments, however, polyploids exhibited fitness advantage, which was conferred by both trait means and adaptive trait plasticities, supporting a ‘jack‐and‐master’ hypothesis for polyploids. Our findings elucidate essential ecological mechanisms underlying polyploid adaptation to heterogeneous environments, and provide an important insight into the prevalence and persistence of polyploid plants.


Fig. S1
Distinct separation between diploid and high-order polyploid Fragaria, with stomatal length as an example.

Common garden soil
Newport beds used a mixture of equal parts beach sand and Bandon fine sandy loam, a low clay content (5-15%) soil derived from the sandy alluvium of coastal marine terraces. Corvallis beds used Chehalis silt clay loam, a richer soil characterized by higher clay content (typically 35-45%). Bend beds used Lundgren ashy sandy loam, a soil characterized by high volcanic ash, glass, and pumice content.

Common garden weather data
We obtained daily data of temperature and rainfall for the three common gardens from different sources: Newport data from the Hatfield Marine Science Center (http://weather.hmsc.oregonstate.edu/weather/weatherproject/archive/), Corvallis and Bend data from AgriMet (https://www.usbr.gov/pn/agrimet/webagdayread.html). We then calculated the monthly mean temperature, monthly rainfall and monthly growing degree days (i.e. the cumulative heat above 10°C) for each garden location, during the course of the field experiment from October 2015 to mid July 2016.

Measurements of leaf functional traits
We collected the largest, fully expanded leaf from each experimental plant in selected beds in each garden. The leaves were scanned using a CanoScan LiDE 220 (Canon, Melville, NY, USA) with an antiglare styrene sheet. We used ImageJ v1.51a (Schneider et al., 2012) to measure leaf area (LA) and the central leaflet width (CLW).
We obtained four leaf punches (each of 6 mm in diameter) from the middle portion of the central leaflet of the collected trifoliate leaf each sample, avoiding the midvein. Two leaf punches were used for measuring stomatal density and stomatal length of the abaxial and adaxial sides; one leaf punch was for measuring SLA; and one was for measuring trichome density and then vein density. When the central leaflet was not large enough for all trait measurements, we obtained two leaf punches from the central leaflet for stomatal density and stomatal length measurements, and one leaf punch from each of the two lateral leaflets for SLA, and trichome density and vein density, respectively.
For SLA estimation, one leaf punch per sample was stored in 96-well microplates (Thermo Fisher Scientific, Hampton, NH, USA), and dried at 65°C for 24 h. Leaf punches were then weighed using a Cahn C-35 microbalance (Thermo Fisher Scientific; with precision of 0.0001 mg). SLA was calculated using the known punch area divided by punch weight.
For stomatal measurements, we used a vinyl polysiloxane impression method to obtain the abaxial and adaxial stomata from leaf punches. First, we mixed the vinyl polysiloxane impression material (Patterson Dental, Pittsburgh, PA, USA) of the base and catalyst, and put the mixture onto a microscope slide. Two punches per sample were placed immediately onto the mixture, one for each side of the leaf. We placed another microscope slide on the top, and pressed slightly and held two slides together using binder clips. After the mixture dried (c. 15 min), the top slide and leaf punches were removed from the mixture using forceps to obtain permanent leaf impression. We applied clear nail polish to the impression and peeled off the impression using clear tapes, and placed it onto a new microscope slide for measuring stomatal density and stomatal length. The abaxial and adaxial stomata were counted using a Leica DM500 microscope (Leica Microsystems, Wetzlar, Germany) under 400× (10 × 40) magnification. Specifically, we counted the total number of stomata within two randomly selected fields of view (FOV) for each side. Stomatal density was calculated as the average number of stomata within a FOV divided by the area of the FOV. We took images of the abaxial and adaxial stomata, and measured the guard cell lengths of up to five stomata of each side and obtained the average stomatal length. As most Fragaria plants only produce stomata on the abaxial side, we only reported abaxial stomatal density and stomatal length.
For trichome density estimation, one leaf punch per sample was stored in 70% ethanol.
We counted the number of trichomes on both sides of a leaf punch under a dissecting microscope. If there were no more than 50 trichomes on one side, we counted all the trichomes. If there were >50 trichomes on one side, we counted trichomes within two randomly selected areas (each area = 1.5 mm × 1.5 mm) of a leaf punch. We summed the abaxial and adaxial trichome density for calculating TD. Leaf punches were then returned to 70% ethanol for subsequent vein density measurement.
Vein density here is defined as the total lengths of minor veins per unit leaf area. We only focused on minor veins, as they account for >80% of the total veins of a leaf and are key to leaf hydraulic capacity and photosynthesis (Sack & Scoffoni, 2013). We followed the protocol of Quantifying Leaf Vein Traits (http://prometheuswiki.org/tikiindex.php?page=Quantifying+leaf+vein+traits), using leaf punches stored in 70% ethanol. We took leaf vein images using a Leica DM500 microscope under 40× magnification, and used ImageJ to record the total lengths of minor veins within a 1 mm × 1 mm area.
The remaining leaf tissue after four leaf punches being taken was dried at 65°C for 48 h, and sent to the Cornell Isotope Laboratory for carbon isotope composition (δ 13 C) and N mass analysis using a Thermo Delta V isotope ratio mass spectrometer and a NC2500 elemental analyzer. Carbon isotope discrimination (Δ 13 C) was calculated using the following formula (Farquhar & Richards, 1984): , ‰, where δ 13 C air equals -8‰.

Plastid phylogeny
The chloroplast nucleotide supermatrix (with 64645 characters), composed of the diploid and polyploid Fragaria taxa in this study (except F. chiloensis ssp. chiloensis) and three other diploid Fragaria taxa, as well as the outgroup Dasiphora fruticosa ssp. floribunda (Fig. S4), was kindly provided by M.S. Dillenberger (Oregon State University). We performed phylogenetic inference using the maximum likelihood (ML) method with the GTR+Γ model in RAxML v.8.0.26 (Stamatakis, 2014). Confidence in node support was determined with 1000 bootstrapping replicates.

Phylogenetic general linear mixed models (PLMMs)
PLMMs were performed to validate our use of nested random effects in LMMs (i.e. populations nested in taxa and taxa in ploidy levels, ploidy/taxon/population; see main text) to control for evolutionary dependence among populations and taxa. Here we used the bifurcating, plastid tree ( Fig. S4) for fitting PLMMs, due to the difficulty of accounting for reticulate evolutionary histories among diploid and polyploid taxa ( Fig. 1) in PLMMs. Owing to the lack of F. chiloensis ssp. chiloensis in the plastid tree, we assumed that it had the same evolutionary history as F. chiloensis ssp. pacifica. We conducted PLMMs using the R package MCMCglmm (Hadfield, 2010), with one functional trait (stomatal length, Fig. S5) and the composite fitness ( To fit MCMCglmm models, we used default priors for predictors (fixed effects) and uninformative priors (V = 1, nu = 0.02) for all random effects and residual variance. Models were run with 200000 total MCMC iterations (burin of 100000, and thinning of 100), and convergence was checked graphically. For Bayesian model comparisons based on the deviance information criterion (DIC), we used the package MuMIn (Bartoń, 2017). Least-squares means of predictors in MCMCglmm models were estimated using the package lsmeans (Lenth, 2016).
For both the functional trait ( Fig. S5) and composite fitness (Fig. S6), LMM Model 2 (the Bayesian version of Model 1) with nested random effects outperformed PLMM Model 4 that only considered phylogenetic relatedness among taxa, but performed as well as PLMM Model 3 that considered both populations nested in taxa and phylogenetic relatedness among taxa.

Fig. S1
Distinct separation between diploid and high-order polyploid Fragaria, with stomatal length (SL) as an example. The least-squares mean of SL and 1 SEM are plotted. Significant differences are only observed between diploids (2n = 2x) and high-order polyploids (2n ≥ 6x). The response variable (SL) was power transformed to improve normality in a general linear mixed model, where the fixed effects included central leaflet width + climatic niche distance + garden + ploidy + ploidy:garden + ploidy:climatic niche distance, and the nested random effects included ploidy/taxon/population. Owing to the distinct separation between diploids and highorder polyploids, and the dominance of the 8x taxa and genotypes (Table S1; also smaller SEM relative to the 6x and 10x here), we defined ploidy level broadly as diploid or polyploid in the main text and all downstream analyses. High-order polyploids Diploids . This worldwide collection of Fragaria was conducted as an international collaborative effort from 2013 to 2014. Each dot represents one population, and the collection data of genotypes within each population are available from the Wild Strawberry website. Briefly, achenes (averagely 70 per plant) were collected from 1-28 plants (mean = 15) of individual populations. For this study, we considered Fragaria that occur in North America, South America, Europe and Japan.

Fig. S3 Climatic niche distances of 72 source Fragaria populations to the common gardens. (a)
The first two principal components of PCA of the 19 bioclimatic variables and elevation estimates of the 72 source Fragaria populations and the three common gardens (the stars). The variables with the largest loadings are indicated by the arrows. (b) The first five PCs, accounting for 94.2% of the variation, were used to calculate Euclidean climatic niche distance between each source population and each garden. This phylogeny reflects only the evolutionary histories of the plastid genome, but not the reticulate histories of the nuclear genome among diploid (2n = 2x) and polyploid (2n ≥ 6x) taxa ( Fig. 1) that are difficult to be incorporated into general linear mixed models for controlling for evolutionary dependence among taxa. This phylogeny included the diploid and polyploid Fragaria in this study (black), and those not (grey). Numbers associated with branches are ML bootstrap support values (%) from 1000 replicates.    Non-parametric Kendall rank correlation coefficient (r) was estimated using the R package psych (Revelle, 2017). Functional trait mean was genotypic trait value averaged across all garden environments.
Non-parametric Kendall rank correlation coefficient (r) was estimated using the R package psych (Revelle, 2017). Functional trait mean was genotypic trait value averaged across all garden environments. Missing r values were due to few data for carbon isotope discrimination and nitrogen content in some taxa.  General linear mixed model (LMM) specification: model <-lmer(Fitness ~ Fixed effects + (1 | Nested random effects) ) Fixed effects: climatic niche distance + garden + ploidy + ploidy:garden + ploidy:climatic niche distance Nested random effects: ploidy/taxon/population The response variable of the LMM was power transformed (power parameter = 0.1). R 2 m , model marginal R 2 representing variance explained by fixed effects; R 2 c , model conditional R 2 representing variance explained by both fixed effects and random effects. General linear mixed model specification: model <-lmer( (Fitness mean) 0.1 ~ Fixed effects + (1 | Nested random effects) ) Fixed effects: climatic niche distance mean + trait plasticity + trait mean + ploidy:trait plasticity + ploidy:trait mean, where climatic niche distance mean and trait mean represent genotypic values averaged across all garden environments Nested random effects: ploidy/taxon/population The response variable (genotypic fitness averaged across all garden environments) of each LMM was power transformed (power parameter = 0.1). R 2 m , model marginal R 2 representing variance explained by fixed effects; R 2 c , model conditional R 2 representing variance explained by both fixed effects and random effects.