Patterns of univariate and multivariate plasticity to elevated carbon dioxide in six European populations of Arabidopsis thaliana

Abstract The impact of elevated carbon dioxide on plants is a growing concern in evolutionary ecology and global change biology. Characterizing patterns of phenotypic integration and multivariate plasticity to elevated carbon dioxide can provide insights into ecological and evolutionary dynamics in future human‐altered environments. Here, we examined univariate and multivariate responses to carbon enrichment in six functional traits among six European accessions of Arabidopsis thaliana. We detected phenotypic plasticity in both univariate and multivariate phenotypes, but did not find significant variation in plasticity (genotype by environment interactions) within or among accessions. Eigenvector, eigenvalue variance, and common principal components analyses showed that elevated carbon dioxide altered patterns of trait covariance, reduced the strength of phenotypic integration, and decreased population‐level differentiation in the multivariate phenotype. Our data suggest that future carbon dioxide conditions may influence evolutionary dynamics in natural populations of A. thaliana.

Characterizing patterns of phenotypic integration can reveal physiological trade-offs that plants may experience in future carbon dioxide conditions (Tonsor & Scheiner, 2007;Ward & Kelly, 2004;Ward & Strain, 1999). For example, carbon supply and acquisition rates can influence water-use efficiency; increased carbon supply may lead to decreased water loss because plants can reduce the density of stomata (i.e., microscopic pores on the surfaces of leaves that regulate gas exchange) without altering carbon acquisition rates (Eamus, 1991;Tonsor & Scheiner, 2007;Ward & Kelly, 2004).
In future carbon dioxide conditions, this may enhance competition of some plants experiencing moderate drought (Woodward, Lake, & Quick, 2002;Lake & Woodward, 2008; although see Perry et al., 2013 for a different perspective). In addition, photosynthesis is typically accelerated in elevated carbon dioxide, resulting in increased carbohydrate production and biomass (Van der Kooij et al., 1999;Makino & Mae, 1999;Teng et al., 2006). However, without correlated increases in nitrogen, elevated carbon supply may reduce photosynthetic rate and generate trade-offs among functional traits. This is due to the correlated acquisition and allocation of carbon and nitrogen, where increased carbon availability leads to increased nitrogen demand (Stitt & Krapp, 1999;Tonsor & Scheiner, 2007;Ward et al., 2000). Therefore, future carbon dioxide conditions may intensify existing trade-offs in nitrogen-deficient environments.
However, evidence also suggests that strong and nonlabile patterns of phenotypic integration may constrain univariate plasticity in plants: tighter integration among traits may reduce the range of variation expressed by individual traits (Gianoli & Palacio-Lopez, 2009;Murren et al., 2015).
An important first step toward predicting short-term responses to elevated carbon dioxide in plants is characterizing univariate and multivariate phenotypic plasticity within and among populations.
Here, we explored population-level differentiation in univariate and multivariate plasticity to elevated carbon dioxide in a set of six functional traits among six European accessions of Arabidopsis thaliana.
We measured total fruit production (reproductive fitness), height (resource allocation to aboveground biomass), rosette diameter (resource allocation during the vegetative phase), biomass (resource allocation during the vegetative and reproductive phases), flowering time (an ecologically important phenological trait), and stomatal density (a functional trait that is crucial for regulating gas exchange).
Six growth chambers (1.2 m length × 0.6 m width × 0.8 m height) were used in this experiment.
Each chamber was fitted with a temperature, carbon dioxide, and humidity sensor (Li-820 and Li-840, LI-COR Biosciences). Carbon dioxide was supplemented using solenoid valve-controlled gas tank regulators (Titan Controls, Inc.) and maintained at 420 ppm in three low carbon dioxide chambers and at 840 ppm, the predicted concentration of atmospheric CO 2 by year 2100 (IPCC, 2014), in three elevated carbon dioxide chambers. Environmental variables were monitored and controlled with Growtronix Software (Tronix Enterprises; growth chamber data logs are available upon request). Plants were grown in 10 × 10 × 8.5 cm pots with standard Pro-Mix soil at 22°C in 16 hr of light per day (150 μmol m −2 s −1 , measured at plant height; HO T12 linear fluorescent tubes were placed approximately 30 cm above plants).

| Experimental design and measurements
To minimize maternal effects, 40 seeds per accession, derived from bulk-propagated populations from the stock center, were randomly assigned to and grown in either the low (420 ppm) or high (840 ppm) carbon dioxide treatment (the maternal population). Six seeds were collected from each maternal plant and grown in the same carbon treatment as the parent (the experimental population). Prior to planting, seeds were imbibed with water and exposed to a five-day cold treatment at 4°C in darkness. In total, the experiment consisted of 2 carbon dioxide treatments × 6 populations × 20 families × 6 replicates = 1,440 plants.
Measurements were taken on the following traits: total fruit production, a measure of reproductive fitness; height, a measure of resource allocation to aboveground biomass; rosette diameter (at bolting), a measure of resource allocation during the vegetative phase; dry weight (after drying at 60°C), a measure of resource allocation during the vegetative and reproductive phases; flowering time, an ecologically important phenological trait; and stomatal density, a functional trait that is crucial for regulating gas exchange and has been shown to respond to changes in carbon availability.
Stomatal density was determined for three fully expanded rosette leaves per plant using the epidermal peel technique described in Ibata, Nagatani, and Mochizuki (2013). Stomata were then visualized using a Nikon Eclipse E200 compound microscope (Nikon Instruments, Inc.) with an integrated digital camera and counted by observing five separate fields of 0.16 mm 2 .

| Data analysis
We used a split plot design where the effect of CO 2 was tested separately from chamber (block)  and p-values were adjusted using a sequential Bonferroni correction when multiple comparisons were carried out. We used three distinct but related multivariate methods to explore variation in phenotypic integration among accessions within and between carbon dioxide treatments. First, we performed principal component analyses and determined trait loadings on leading eigenvectors. Second, we estimated eigenvalue variance to assess the strength of phenotypic integration (Armbruster, Pélabon, Bolstad, & Hansen, 2014;Pavlicev, Cheverud, & Wagner, 2009;Wagner, 1984). Stronger integration results in the concentration of variation in leading eigenvalues, which in turn leads to higher eigenvalue variance. In contrast, weaker integration produces similar eigenvalues with lower eigenvalue variance.
Third, we performed pairwise comparisons of covariance matrices between unique population-treatment groups using common principal components analysis (CPCA; Phillips & Arnold, 1999). CPCA is a multivariate statistical method used to detect structural similarities among covariance matrices. It is based on Flury's (1988) hierarchy and has been adapted to quantitative genetic data by Phillips and Arnold (1999). In contrast to approaches that test for only equality, CPCA tests a series of nested models of structural similarity between covariance matrices. The model hierarchy begins with the assumption of unrelated covariance structure. This is followed by a series of increasingly constrained models that include partial CPCs, full CPCs, proportionality, and equality. This approach enables the detection of subtle differences in covariance structure. Using the "jump-up" approach of CPCA (Phillips & Arnold, 1999), we indicated the best-fitting model for each comparison as the one below the point in the hierarchy where CPCA found statistically significant differences between matrices (i.e., the best model in the hierarchy is the one below the rejected model).

| RE SULTS
Analysis of variance revealed significant population-level variation in all traits and significant within-population, family-level variation in only one trait (fruit number; Table 1). Elevated carbon dioxide had significant effects on trait means (Table 1). Stomatal density decreased by 11.28%, and the other five traits increased by 10.87%-43.61% (Table 2). Only one trait (flowering time; Table 1) showed significant variation in plasticity among accessions. In contrast, we did not detect significant family-level variation in plasticity (i.e., GxE interaction within populations) in any trait (Table 1).
Correlation analysis showed 14 out of 15 significant correlations in the 420 ppm treatment and 12 out of 15 significant correlations in the 840 ppm treatment (Table 3). We found three significant negative correlations in 420 ppm treatment and one significant negative correlation in 840 ppm treatment. The correlation between rosette diameter and stomatal density changed from weak negative (r 2 = −0.08; p < 0.05) in 420 ppm to moderately strong positive (r 2 = 0.49; p < 0.05) in 840 ppm (Table 3). The correlation between weight and flowering time changed from positive (r 2 = 0.22; p < 0.05) in 420 ppm to negative (r 2 = −0.34; p < 0.05) in 840 ppm (Table 3).
Elevated carbon dioxide resulted in decreased eigenvalue variance, a measure of the strength of phenotypic integration, in all accessions except Sidmouth, UK (Figure 1). Principal components analyses revealed that the first two components accounted for 62.63% and 57.24% of the variance in low and high carbon treatments, respectively (Table 4). Trait loadings showed that in 420 ppm treatment, principal component 1 (PC1) was influenced mainly by variation in weight, rosette diameter, height, and stomatal density, and PC2 was influenced mainly by flowering time and fruit number (Table 4). In contrast, in 840 ppm treatment, fruit number and weight loaded more heavily onto PC1 and PC2, respectively (Table 4)

| D ISCUSS I ON
The unprecedented rise of atmospheric carbon dioxide-a key driver of global climate change-represents a significant ecological

| Univariate plasticity to elevated carbon dioxide
We detected significant population-level variation in all traits and significant family-level variation in only fruit production. These patterns are consistent with previously documented quantitative variation in A. thaliana (Koornneef, Alonso-Blanco, & Vreugdenhil, 2004;Kuittinen, Mattila, & Savolainen, 1997;Pigliucci, 2002), a mostly self-fertilizing species that is expected to harbor more genetic variation among than within populations. In addition, while elevated carbon dioxide-induced plastic responses in all traits, we did not detect significant variation in plasticity (i.e., genotype by environment interaction, or G × E) in any trait except flowering time. Plasticity to TA B L E 3 Correlations among six traits of Arabidopsis thaliana in two carbon dioxide treatments (420 ppm = below diagonal; 840 ppm = above diagonal)  Springer & Ward, 2007).
Stomatal density is sensitive to a wide range of CO 2 concentrations (including preindustrial levels), and responses in this trait have been shown to vary in a genotype-dependent manner in Arabidopsis (Lake & Woodward, 2008;Ward & Kelly, 2004;Ward & Strain, 1999;Woodward et al., 2002). We found an approximate 11% decrease in stomatal density in elevated carbon dioxide, which is consistent with some previous findings for this trait under carbon enrichment (Lake & Woodward, 2008;Woodward et al., 2002). This result supports the notion that plants in carbon-rich environments can afford to reduce stomatal densities while maintaining carbon intake and water-use efficiency (Eamus, 1991;Tonsor & Scheiner, 2007). Our experiment was not designed to measure photosynthetic rates or water-use efficiency, but future empirical work may consider quantifying these traits along with stomatal densities in a multivariate context. We suggest that the patterns of integration among these traits may likely be altered by elevated carbon dioxide, and trade-offs among photosynthetic rate, biomass, and stomatal densities may arise, especially in nitrogen-deficient environments due to shared resource allocation (see Introduction ;Stitt & Krapp, 1999).
In addition to significant variation in flowering time responses described above, based on previous work we expected to find significant genetic variation in plasticity of fitness (Ward & Strain, 1999), stomatal density (Woodward et al., 2002), and limited genetic variation in the plasticity of traits associated with the production of biomass (rosette diameter, height, and biomass; Ward & Strain, 1999;Ward & Kelly, 2004). However, our results are, overall, in line with work showing significant plastic responses to elevated CO 2 in morphological, physiological, phenological, and fitness-related traits (reviewed in Leakey & Lau, 2012;Ward & Kelly, 2004;Ward & Strain, 1999;Wieneke, Prati, Brandl, Stöcklin, & Auge, 2004), but little or no genetic variation in phenotypic plasticity to elevated CO 2 (Lau, Shaw, Reich, Shaw, & Tiffin, 2007;Leakey & Lau, 2012;Roumet, Laurent, Canivenc, & Roy, 2002;Volk & Körner, 2001; although see Lindroth, Roth, & Nordheim, 2001). Specifically, among studies focused on variation in plasticity to elevated CO 2 in biomass and fitness (seed number, fruit number, seed weight, or relative growth rate), less than a third have found statistically significant GxE in these traits: Lau et al. (2007) and Roumet et al. (2002) report 7/21 and 11/39 studies, respectively, that have found significant GxE in such traits.
Furthermore, shifts toward weaker phenotypic integration are predicted to facilitate independent variation in traits and promote evolutionary potential in univariate dimensions. In turn, changes in phenotypic integration may alter patterns of selection and evolutionary outcomes in A. thaliana. Evidence suggests that elevated carbon dioxide can, in fact, relax competition-mediated selection (Lau et al., 2014) and influence the strength of selection on carbon assimilation in Arabidopsis (Tonsor & Scheiner, 2007).
We found some degree of population-level differentiation in phenotypic integration in low carbon dioxide treatment in both the UK and PO regions, as indicated by pairwise CPCAs of trait covariance matrices (Figure 2). It is well known that A. thaliana is a highly selfing species and ecotypes are differentiated in their life history characteristics (Koornneef et al., 2004;Pigliucci, 2002 In contrast to low carbon dioxide treatment, we did not detect differentiation among trait covariance matrices among accessions in elevated carbon treatment in either the UK or PO regions ( Figure 2). Although we used a unique set of statistical approaches (correlation analysis, eigenvalue variance analysis, PCA, CPCA), our finding that trait integration is sensitive to changes in carbon availability is in general agreement with previous experiments that used principal components analyses (Bidart-Bouzat et al., 2004) and structural equation modeling (Tonsor & Scheiner, 2007) to characterize the effects of carbon enrichment on phenotypic integration in A. thaliana. Therefore, this finding appears to be robust to both variation in traits measured and differences among multivariate techniques used to analyze patterns and magnitudes of integration.
In addition, our predictions that elevated carbon dioxide will likely mitigate trade-offs and relax constraints on adaptation are in line with a previous study that also found significant shifts in patterns of integration, as well as partial decoupling between traits (carbon assimilation and transpiration rates) and weaker selection on carbon assimilation in A. thaliana accessions grown across a CO 2 gradient (Tonsor & Scheiner, 2007). Finally, our findings are overall consistent with the broader idea, supported by a large body of work, that rapid changes in the environment can have significant impacts on phenotypic integration (Bidart-Bouzat et al., 2004;Chevin et al., 2010;Lind et al., 2015;Marroig & Cheverud, 2001;Moczek et al., 2011;Murren, 2002Murren, , 2012Pigliucci & Preston, 2004;Plaistow & Collin, 2014;Price et al., 2003;Schlichting, 1989;Sgrò & Hoffmann, 2004;Wund, 2012 Future studies that combine ecological complexity with populationlevel multivariate analyses of phenotypic plasticity and integration could deepen our understanding of evolutionary potential of plants in the face of climate change.

ACK N OWLED G M ENTS
Discussions with and edits by Mike Bell, R. Geeta, and Massimo Pigliucci greatly improved the final version of this manuscript. We thank Cathy D'Agostino, Natalie Odynocki, Lily Sarrafha, and Tianyu Shang for help with empirical work. We also thank two anonymous reviewers for their thoughtful feedback on a previous version of the manuscript.

CO N FLI C T O F I NTE R E S T
None declared.

AUTH O R CO NTR I B UTI O N S
M.J. designed the study, collected data, analyzed data, and wrote the manuscript. B.C. collected data and analyzed data.

DATA ACCE SS I B I LIT Y
Data are available on Dryad: https ://doi.org/10.5061/dryad. c0k235b.