Plant adaptation or acclimation to rising CO2? Insight from first multigenerational RNA‐Seq transcriptome

Atmospheric carbon dioxide (CO2) directly determines the rate of plant photosynthesis and indirectly effects plant productivity and fitness and may therefore act as a selective pressure driving evolution, but evidence to support this contention is sparse. Using Plantago lanceolata L. seed collected from a naturally high CO2 spring and adjacent ambient CO2 control site, we investigated multigenerational response to future, elevated atmospheric CO2. Plants were grown in either ambient or elevated CO2 (700 μmol mol−1), enabling for the first time, characterization of the functional and population genomics of plant acclimation and adaptation to elevated CO2. This revealed that spring and control plants differed significantly in phenotypic plasticity for traits underpinning fitness including above‐ground biomass, leaf size, epidermal cell size and number and stomatal density and index. Gene expression responses to elevated CO2 (acclimation) were modest [33–131 genes differentially expressed (DE)], whilst those between control and spring plants (adaptation) were considerably larger (689–853 DE genes). In contrast, population genomic analysis showed that genetic differentiation between spring and control plants was close to zero, with no fixed differences, suggesting that plants are adapted to their native CO2 environment at the level of gene expression. An unusual phenotype of increased stomatal index in spring but not control plants in elevated CO2 correlated with altered expression of stomatal patterning genes between spring and control plants for three loci (YODA, CDKB1;1 and SCRM2) and between ambient and elevated CO2 for four loci (ER, YODA, MYB88 and BCA1). We propose that the two positive regulators of stomatal number (SCRM2) and CDKB1;1 when upregulated act as key controllers of stomatal adaptation to elevated CO2. Combined with significant transcriptome reprogramming of photosynthetic and dark respiration and enhanced growth in spring plants, we have identified the potential basis of plant adaptation to high CO2 likely to occur over coming decades.


Introduction
Since industrialization, global atmospheric CO 2 concentrations ([CO 2 ]) have risen steadily from ca. 280 ppm to present-day concentrations of ca. 400 ppm and are predicted to rise to between 700 and 900 ppm. Although considerable evidence exists to suggest that CO 2 has acted as a driver for plant evolution over past geological time (Parmesan & Yohe, 2003), evidence of a role for future [CO 2 ] in determining plant adaption is limited (Ward & Kelly, 2004). Until recently, the small number of studies provided inconclusive findings (Leakey & Lau, 2012). This is surprising given the large impacts of elevated [CO 2 ] on plant traits including those linked to productivity (Ainsworth & Long, 2005) and fitness (Springer & Ward, 2007). On the other hand, phenotypic plasticity has already been observed in plant response to [CO 2 ], for example in a QTL mapping population, where genotypic difference in response to elevated [CO 2 ] resulted in several G 9 E interactions and response QTL identified in Populus trichocarpa T. & G 9 P. deltoides Marsh (Rae et al., 2007). Multigenerational exposure to elevated [CO 2 ] in the model photosynthesizing organism Chlamydomonas (Collins & Bell, 2006) has also revealed subtle changes in photosynthetic capacity and cell growth characteristics from a high CO 2 spring, but no underpinning genetic and genomic data were available to understand these responses further. Very recently, Grossman & Rice (2014) and Horgan-Kobelski et al. (2016) also identified reaction norm differences when contrasting populations were exposed to ambient and elevated [CO 2 ] including Bromus madritensis exposed for several years to elevated [CO 2 ] in a free-air CO 2 enrichment (FACE) experiment and then subjected to either ambient or elevated CO 2 in a controlled environment. These experiments suggest strongly that plants are likely to show adaptive responses to this important facet of the changing climate but that as yet, these responses are not fully understood.
Changes may be subtle and therefore require appropriate experimental systems for further investigation, which have not to date been forthcoming. Leakey & Lau (2012) suggest that one of the limitations to our current understanding is the necessary emphasis on simplistic and short-term experiments. Few multigenerational exposures to elevated CO 2 have been undertaken. Coupled with this is the limited genomic and genetic characterization of plants from long-term experiments. The availability of new molecular tools, in particular high-throughput and inexpensive RNA and DNA sequencing, suggests that this bottleneck may now be addressed using previously impossible approaches that combine phenotyping, functional genomic and population genetic analysis in nonmodel systems.
Natural elevated CO 2 springs occur where CO 2 is emitted from vents, increasing local [CO 2 ] to more than double that only a few metres away (Bettarini et al., 1999). It is predicted that [CO 2 ] in such springs has been elevated for hundreds to thousands of years (Bettarini et al., 1999;Cook et al., 2008), suggesting that plants that live within the springs have the potential to show adaptation to these conditions, as opposed to being acclimated in the short-term and expressing phenotypic plasticity. To date, no molecular genomic and genetic analyses have been undertaken on these putatively adapted plants, although they provide a valuable resource for improved understanding.
When plants are exposed to elevated [CO 2 ], gene expression is altered (Leakey et al., 2009a;Tallis et al., 2010) and this indicates that changes in gene regulation could be a prominent mechanism underpinning adaptation to elevated [CO 2 ]. Past transcriptomic studies have identified photosynthesis (Taylor et al., 2005;De Souza et al., 2008), respiration (Leakey et al., 2009a) and leaf development (Ainsworth et al., 2006) as plant processes where gene expression is sensitive to elevated CO 2 . Thus, a link has been made between plant traits responsive to elevated CO 2 and changes in the expression of key genes including genes coding the small subunit of RuBisCo and cell wall loosening enzymes such as XTH, linked to plant growth. Such gene expression changes, linked to important adaptive traits, represent phenotypic plasticity and could offer a clue to the targets of selection during long-term (multigenerational) adaptation. It is also possible that in addition to gene regulation, mutation acts to give rise to novel locally adapted alleles. To investigate these possibilities, we collected seed of Plantago lanceolata L. from a population located in a CO 2 spring where plants have been subjected to elevated [CO 2 ] over multiple generations, the spring population (S) and from a nearby control population (C) site with ambient [CO 2 ]. Progeny were subjected to either ambient or elevated [CO 2 ] in controlled conditions, followed by detailed phenotypic and transcriptomic (RNA-Seq) characterization, where traits were chosen that underpinned growth and leaf development as these are known to be linked to plant fitness and have previously been shown to be responsive to atmospheric [CO 2 ]. Much research has been carried out at these and other springs in the past, on photosynthesis, growth, stomatal development and leaf chemistry, but no functional or population genomic studies have previously been reported and only limited population genetic data are available (Nakamura et al., 2011). From the RNA-Seq data † (gene expression and DNA polymorphism), it was possible to quantify (1) differences in gene expression due to acclimation to elevated [CO 2 ], (2) differences in gene expression likely due to evolutionary adaptation to the [CO 2 ] spring, (3) the strength of genetic differentiation between the spring and control plants and finally (4) the link between these data to phenotypes from spring and control plants exposed to either ambient or elevated [CO 2 ].

Materials and methods
Plant material and sampling site Plantago lanceolata L. (Plantaginaceae) seeds were collected from a natural CO 2 spring in Bossoleto, Italy (Lat. 43°17 0 , Long. 11°35 0 ), where the atmospheric [CO 2 ] is maintained during the day at an average concentration between 600 and 1200 lmol mol À1 whilst unstable nocturnal conditions are often >7500 lmol mol À1 (Bettarini et al., 1999). Plantago lanceolata growing in the elevated [CO 2 ] environment of the spring (spring site, S) were compared with the same species growing outside the spring in ambient [CO 2 ] (control site, C, 200 m from spring). Seeds of P. lanceolata were collected from nine randomly selected maternal plants in both sites on 12 May 2008 and stored at 5°C P. lanceolata is a short-lived perennial rosette herb that requires long-day conditions to induce flowering. It is self-incompatible and therefore an obligate outbreeder. It has been suggested that P. lanceolata is capable of rapid evolutionary change, following exposure to contrasting environmental challenges. Controlled environment exposed to elevated [CO 2 ] Seeds from each family were sown in a glasshouse at The University of Southampton, Boldrewood campus. In order to reduce maternal effects when assessing population response to elevated [CO 2 ], ten families (sibs from one maternal plant) from each site were established, isolated in a muslin tent and allowed to set seed. The way that this material was used for the phenotypic measurements and gene expression data is detailed in Table S1. Seeds from each family were sown in 20-cm-diameter pots in compost (John Innes potting compost No. 2; John Innes Manufacturers Association, Norwich, UK) and topped with 1-cm fine sand. Three weeks after germination, plants were moved into one of eight CO 2 chambers, with one member of a pair of siblings from each family allocated to ambient (A) and one to elevated (E) treatments (Fig. 2a). Four A chambers were maintained at a [CO 2 ] of 410.63 AE 33.74 ppm, and four E chambers maintained at e [CO 2 ] of 718 AE 46.81 ppm (flow rate 3.4 m s À1 and PAR at 104-134 lmol m À2 s À1 ). The position of the pots was randomized within and between each chamber every 2 days, and treatments were swapped between chambers weekly to remove chamber effects (Warwick & Taylor, 1995).
In this study, traits were selected for analysis that are known from previous research to be sensitive to elevated atmospheric CO 2 , in particular biomass production and also leaf size and the determinants of leaf sizecell size and number. In addition, stomatal development was also assessed from stomatal density (SD) and stomatal index (SI) measurements because both are known to be sensitive to elevated CO 2 (Engineer et al., 2014). Although photosynthesis, stomatal conductance and respiration are also sensitive to elevated CO 2 , these are not reported here, although our current research is quantifying family response of P. lanceolata to elevated [CO 2 ] (K. Staniland, J. Saban, M.A Chapman, G. Taylor, personal communications).

Biomass traits
Whole plants were harvested on day 124 of the experiment. The second or third young expanding leaf was harvested into liquid nitrogen and stored at À80°C. A mature leaf of each plant was scanned (HP Scanjet 4850; Hewlett-Packard Development Company, Palo Alto, CA, USA) and measured using ImageJ (Image J 1.42q; Wayne Rasband, Bethesda, MA, USA) for single leaf area. The whole plant was cut at 1 cm above ground, and all leaves were laid on a white scaled A3 paper for imaging. Mature leaves were weighed separately after drying at 75°C for 2 days or until leaves no longer changed weight.

Stomatal traits
Cell imprints were taken from the abaxial surface of a mature leaf from each plant at final harvest. Nail polish was painted onto the middle of the abaxial surface of the leaf and imprints were taken as described previously (Potvin & Tousignant, 1996). The number of stomatal and epidermal cells was counted on images under fields of view of 1.6 mm 2 at 109 magnification. Epidermal cell size (ECS) was measured using twenty randomly selected cells and guard cells in each image. Specific leaf area, SI, SD and epidermal cell number per leaf were calculated as described previously (Taylor et al., 2003).

Statistical model
Morphology measurements were tested by generalized linear models in SPSS 16 (SPSS Inc, Chicago, IL, USA) with the model: Response = Location + Treatment + Treatment | Location + Treatment | Family (Location) + Location | Chamber (Treatment). Traits were log-transformed before analysis if normalization tests were rejected (P < 0.05) with a Kolmogorov-Smirnov test and planned post hoc Sidak's test was applied on the traits which were significantly influenced by location and treatment interaction.

RNA-Seq analysis of the transcriptome
For transcriptomic analysis, the expressed portion of the genome of 24 samples [six each of spring ambient (SA); spring elevated (SE); control ambient (CA); control elevated (CE); Fig. 1b] were sequenced and gene expression confirmed on a subset of genes using real-time PCR (Fig. S1). RNA was extracted from young leaves using a CTAB-based protocol (Taylor et al., 2005) and quantified using a Bioanalyser (Agilent Technologies, Santa Clara, CA, USA). Library preparation with poly(A) mRNA selection was carried by oligo d(T) magnetic beads on each total RNA sample. mRNA was fragmented at 94°C using divalent cations, and each fragment of mRNA was synthesized to cDNA. After double-strand cDNA reparation, poly(A) was added to the 3 0 ends of each fragment. Sequencing adapters were then added to each fragments and amplified by PCR. Sequencing took place at the Instituto di Genomica Apllicata using an Illumina HiSeq 2000 (Illumina Inc., San Diego, CA, USA), producing a 100-bp paired-end reads. Following basic quality control, the Trinity pipeline (Grabherr et al., 2011) on the Iridis 4 supercomputer at the University of Southampton Supercomputing facility was used for de novo transcriptome assembly and downstream analyses as follows. Further details and Perl scripts are available from the authors upon request. First, to reduce the overall number of reads that were used in the initial assembly, reads from eight libraries (two from each population/treatment) were normalized using normalize_by_kmer_coverage.pl with maximum coverage set to 30. The resultant libraries were assembled into a de novo transcriptome using Trinity.pl with min_kmer_cov set to 2 to reduce the proportion of error-containing kmers, and taking into account the orientation of paired reads. This resulted in a transcriptome comprising 154 179 transcripts from 43 908 components. Transcripts represent different forms of the same component, and components are loosely comparable to genes; however, two components may represent different portions of the same gene.
Individual libraries (non-normalized) were mapped back to the transcriptome using run_RSEM_align_n_estimate.pl and then converted to .bam files for sequence analysis. Because of the high number of transcripts, most of which had very low expression, we carried out differentially expressed (DE) analysis on the sum of reads that mapped to each component (gene). Counts were converted to TMM-normalized FPKM using abundance_estimates_to_matrix.pl in Trinity, that is to standardize by gene length and by total number of reads per sample. A principal component analysis (PCA) was performed on the expression data using the prcomp function of R (version Rx64 3.0.1; R Core Team, Vienna, Austria). Data used for PCA were TMM-normalized FPKM and were log 10 -transformed. In order to interrogate gene expression data to full effect and in a robust manner, two separate analyses were completed for DE across the four treatments. edgeR was used to identify DE genes, given its value in using Poisson rather than normally distributed data that may be beneficial for counts such as in RNA-Seq data sets (Robinson et al., 2010). It was used in run_DE_analysis.pl (from within Trinity) and DE determined between the four location/treatments using ana-lyze_diff_expr.pl and a 5% FDR cut-off. The second analysis using t-tests and 5% FDR (using Bonferroni correction) was also then completed to determine whether the expression values differed significantly across treatments or source populations. In a comparison of these two approaches, the trends of DE were similar in both analyses and in general more than 70% of DE data with the same genes were significant in both analyses. Presented here are t-tests with FDR at 5%. To avoid missing biologically important differences, two-way ANOVA with FDR at 5% are given for subgroups of the data to show significantly different genes within subgroups.
For further detailed functional genomic analysis, DE gene data were explored using MAPMAN (Thimm et al., 2004) and AGRIGO (Du et al., 2010). Reads were first assembled using de (a) Growth and morphological responses of control and spring plants exposed to either ambient or elevated CO 2 . (a) Reaction norms for above-ground biomass, (b) above-ground biomass, (c) epidermal cell size, (d) single leaf area, epidermal cell number per leaf, (e) stomatal index (f) specific leaf area. Mean (AEstandard error) and the ANOVA results for location (S/C), treatment (A/E) and the interaction term are shown with complete statistical analysis given in ST1. Stomatal imprints from either control (g) or spring (h) site, grown at e[CO 2 ] (700 lmol mol À1 ). (i) Percentage change in biomass, cell and stomatal traits of plant grown in either a[CO 2 ] (~400 lmolÁmol À1 )or e[CO 2 ] (700 lmolÁmol À1 ) (percentage change = (elevated-ambient)/ambient)*100).
novo assembly (Grabherr et al., 2011), and then, the resulting transcripts were mapped onto the best Arabidopsis thaliana (AT) genes using BLASTN (E < 10 À5 ). If multiple Plantago transcript matches were made to a single AT gene, then average expression was used for analysis to ensure all isoforms were represented. The MAPMAN analysis used DE genes identified from ANOVA (FDR < 0.05). For AGRIGO, a Fisher significance test at P < 0.05 was used for DE genes and analysis restricted to GO categories with at least five mapping entities, enabling over-enriched functional categories to be identified (Figs S3, S4). AT genes known to be involved in the regulation of stomatal numbers and patterning (ER, ERL1, YODA, MYB 88, FAMA, SPCH, CDKB1;1, SCRM2, EPF1, BCA1, BCA4) were identified from the literature (Torii, 2015) and mapped to transcripts as described above. For ER, BCA1 and the RuBisCo small subunit, DE data were confirmed using real-time PCR (RT-PCR, Fig. S1, Table S3), as described in Supporting information.

Population genomics
Because of the relatedness between pairs of plants subjected to A and E [CO 2 ], this analysis only considered six plants from different mothers from each of the C and S sites. Transcript sequence sorted BAM files from Trinity were parsed through SAMtools (Li et al., 2009) to (1) align the reads to the reference de novo assembly, (2) create consensus (with heterozygous bases encoded with IUPAC code) for each locus and (3) to save sequences in fasta format, using the following settings: minimum base quality 13, minimum mapping quality 2 and minimum read depth 3. Resultant fasta files were imported into PROSEQ (Filatov, 2009) and aligned to the reference de novo assembly. Heterozygous bases were phased by exporting the alignments, running each through FASTPHASE (Scheet & Stephens, 2006) and importing the resultant phased files back into PROSEQ. Perl scripts are available on request. The following analyses only considered the 12 428 loci with sequences from all 12 individuals present. Nucleotide diversity, p, and Tajima's D (Nei & Li, 1979), a measure of neutrality, were calculated in PROSEQ (Filatov, 2009) and compared between species using Mann-Whitney U-tests in PAST (http://folk.uio.no/ ohammer/past/). Per locus differentiation between the C and S populations was estimated using Wright's FST (Wright, 1949) within PROSEQ. PROSEQ was also used to output files for use with BAYESCAN (Foll & Gaggiotti, 2008), a program designed to identify loci with high interpopulation differentiation. The program was run with 150 000 iterations, discarding the first 50 000 as burn-in. Outlier loci based on a 5% FDR were identified.

Growth and morphological changes in response to elevated [CO 2 ]
The morphological data present clear evidence of differences between the spring and control plant responses to elevated [CO 2 ] (Fig. 1a-f). For aboveground biomass (Fig. 1a, Table S2), a consistent treatment effect was apparent with larger plants in elevated [CO 2 ] for both spring and control plants. A significant Location 9 Treatment effect (P < 0.05) was noted for both SI and single leaf area (Table S2). Stomatal index was 5.21% greater in spring plants (Fig. 1g, h), but showed no difference for control plants exposed to either ambient or elevated [CO 2 ]. Single leaf area was 5.51% smaller in spring plants and 12.23% greater in control plants when grown in elevated [CO 2 ] and compared to ambient [CO 2 ]. Epidermal cell number per leaf exhibited a significant treatment response (P < 0.05), with an increase for control plants in response to elevated [CO 2 ]. Single leaf dry mass showed a significant location effect (P < 0.05), with the spring plants having consistently greater single leaf dry mass across CO 2 treatments (average 20% greater). These data are plotted as norms of reaction, revealing, overall, considerable evidence that phenotypic response to elevated [CO 2 ] varied between the spring and control populations. This is apparent from the percentage response to elevated CO 2 data (Fig. 1i).

Gene expression changes in response to elevated [CO 2 ]
Although differential gene expression was apparent for both spring and control plants in response to elevated [CO 2 ] (Fig. 2b- (Fig. 2d). This response was further characterized by the PCA, where PC1 clearly separates CA and CE plants, but not SA and SE plants (Fig. 2c). This figure suggests that the transcriptomic response to elevated [CO 2 ] for the control plants was much larger than that of the spring plants and that for spring plants, gene expression appears 'adapted' to the elevated CO 2 environment of the elevated CO 2 treatment. Gene expression for specific metabolic pathways was explored using MAPMAN (Fig. 3), which gives a critical insight into complex plant processes such as photosynthesis and how groups of genes that control steps in many different pathways may be effected. MAPMAN Fig. 3a). Over 80% of the genes related to photosynthesis were significantly DE (FDR < 0.05) using two-way ANOVA (Fig. 3a, defined fully in Table S3). For example, RuBisCo small subunit gene expression was overall 28.6% lower for control plants in elevated [CO 2 ] compared with ambient [CO 2 ] [significant using two-way ANOVA and FDR < 0.05 (Fig. 3, Table S3)]. Spring plants, in contrast, showed only minimal difference in gene expressionfor example, a 5.3% reduction in expression in response to elevated [CO 2 ] across the same RuBisCo small subunit transcripts. Respiration followed a similar pattern, with control plants showing on average, a 119% reduction in expression across the 'TCA cycle and mitochondrial electron transport' category in response to elevated [CO 2 ], whilst for spring plants expression was 7.4% increased for this MAPMAN category in response elevated [CO 2 ] (Fig. 3a, b). AGRIGO analysis (Figs S3 and S4) supported these findings, revealing photosynthetic and respiration-related functional GO categories as significantly overrepresented in the DE list (FDR < 0.05). Taken together, these MAPMAN  The stomatal patterning pathway (Fig. 4) is of special interest due to the contrasting stomatal patterning phenotypes of spring and control plants exposed to ambient and elevated [CO 2 ] (Fig. 1)  patterning has already been elucidated. In our data set, we identified thirteen homologues of genes previously documented as related to stomatal patterning, seven of which were DE and one (FAMA) DE at the 10% level of probability for both treatment and location (using an FD < 0.05; Fig. 4 where DE homologues are shown). The comparison of plants grown in ambient and elevated [CO 2 ] -ANOVA treatment effect, revealed four DE stomatal patterning genes (FDR < 0.05), with tentative evidence for DE (FDR < 0.1) for a further two (Fig. 4). The beta carbonic anhydrase (BCA) double mutant bca1 bca4 has a 20% increase in stomatal numbers in elevated [CO 2 ] (Hu et al., 2010), and here for both spring and control plants, a putative orthologue of BCA1 was expressed significantly less in elevated [CO 2 ], suggesting a role in determining stomatal patterning in Plantago in elevated [CO 2 ]. Of greater biological significance however is the comparison of two genes (SCRM2 and CDKB1;1) that are expressed during the asymmetric divisions early in the stomatal patterning pathway (Torii, 2015), both of which were expressed significantly more in elevated CO 2 in spring plants relative to ambient-grown spring plants.

Population genetics of the spring and control plants
Population genomic analysis revealed essentially no genetic differentiation between the populations, with mean F ST across all loci being only 0.045 AE 0.064. Although we acknowledge a relatively small sample size (12 alleles per population), this is unlikely to bias the overall F ST downwards. A test for outlier loci (i.e. those with greater than expected genetic differentiation between the two populations, see methods) failed to uncover any such loci. Whilst small sample size could again affect our ability to detect outlier loci, a comparison of outlier detection methods, Narum & Hess (2011), showed Bayescan to have the lowest type II error (i.e. the smallest chance of not identifying an actual outlier). Sequence polymorphism (p) was significantly higher in the control population than in the spring population (p Control = 0.0127 AE 0.0109 (SD) vs. p Spring = 0.0103 AE 0.0099; paired t-test; t = 37.28; P < 0.0001), and Tajima's D was also significantly different between the two populations (D Control = À1.2476 AE 0.8718 [SD] vs. D Spring = À0.6074 AE 1.0151; paired t-test; t = 73.45 P < 0.0001).

Discussion
Understanding plant adaptation to a rapidly changing climate remains a high scientific priority, given the implications for food security and in the conservation and management of a wide range of ecosystems (Mawdsley et al., 2009;Branca et al., 2013). Despite this, only limited research is available on likely multigenerational responses to increasing atmospheric [CO 2 ]. This was confirmed in a recent meta-analysis to assess plant fitness in relation to climate change (Anderson, 2016), where a very limited data set was identified that described studies on elevated CO 2 fitness traits, in contrast to those from increased temperature experimental treatments. Thus, the research in this present study addresses an important research question, as the likelihood that atmospheric CO 2 can act as a selective pressure driving evolution seems high, given the direct link with photosynthesis, carbon balance and indirect link to other fitness-related traits including growth and reproduction.

Changes in growth and morphology of spring plants may indicate pre-adaptation to elevated [CO 2 ]
Growth and morphological changes have been observed previously when plants are grown in elevated [CO 2 ], including increased plant biomass, leaf size and cell size and number (Taylor et al., 1994;Ainsworth and Long, 2005). Here, we found that spring and control plants behaved differently in response to elevated [CO 2 ], indicative of local adaptation, in particular the reaction norms differed significantly for leaf area, SI and leaf morphology traits. Control plants showed increased leaf area, epidermal cell number and total biomass in response to elevated [CO 2 ], in agreement with previous research. In contrast, spring plants often showed a different response to the elevated [CO 2 ] treatment (Fig. 1), although the interaction term was often not significant, despite these contrasting responses. This likely reflects the power of the experimental design where 9-10 contrasting families were used in the analysis representing genetic and phenotypic diversity. Despite this, the data here strongly suggest altered phenotypic plasticity between spring and control plants, as observed recently for stomatal conductance in a long-term adapted FACE population of Bromus (Grossman & Rice, 2014). Although past studies show that the majority of plants exposed to elevated [CO 2 ] have reduced numbers of stomata (measured as SI or density (Hetherington & Woodward, 2003), here we found SI was significantly increased for spring plants in response to exposure to elevated [CO 2 ] (Fig. 1h, i), a finding reported previously in P. lanceolata (Marchi et al., 2004) Up to 30% of species within certain subclasses show increased and not decreased SD and index in response to elevated [CO 2 ] (Woodward & Kelly, 1995), as well as cases of greater SD within naturally high CO 2 springs, and that the finding here for P. lanceolata was observed in three separate experiments (data not shown) and is thus robust and should be considered further. The consequences of increased stomatal numbers for adaptation of leaf gas exchange are intriguing as a greater SI in this study was often associated with stomatal 'clumping', where several stomatal pores were found together that did not comply with the 'one cell' rule, where pores are consistently separated by at least one cell (Dow & Bergmann, 2014). Stomatal clumping leads to reduced stomatal conductance as gaseous concentration shells are joined, which may reduce flux by 5-15% (Lehmann & Or, 2015). On the other hand, a greater number of wellspaced stomata may enhance leaf stomatal conductance and water loss and so it remains a point of speculation as to overall impact of this phenotype on leaf gas exchange and plant adaptation. Further research is ongoing to elucidate these functional gas exchange traits and how they link to underpinning stomatal guard cell gene expression.

Few gene expression changes in spring plants compared with control plants when exposed to elevated [CO 2 ]
Our data suggest for spring plants gene expression changes were modest in response to elevated CO 2 , with somewhat larger differences found for control plants exposed to the different CO 2 treatments, illustrated in the PCA where CA and CE gene expression data sets were separated but those for SA and SE were not. The smaller number of DE loci and lack of separation in gene expression data sets for the spring plants could result from adaptation of the spring plants to elevated [CO 2 ], in line with a number of limited gross phenotypic changes observed for the spring plants (Fig. 1). Overall, gene expression changes in response to [CO 2 ] were moderate when compared to the large number of differences observed between control and spring plants irrespective of CO 2 treatment (Fig. 2C). Plants exposed to elevated [CO 2 ] often respond by downregulation of the photosynthetic machinery (Ainsworth & Long, 2005), altered respiration (Leakey et al., 2009b) and an associated increase in plant growth and biomass production (Ainsworth & Long, 2005). For control plants, we observed these responses, with lower RuBisCo gene expression, despite a 20% greater biomass at elevated CO 2 (Figs 1, 3, 5 and Table S2). In addition to this, several genes coding for the respiratory machinery of the cell exhibited reduced expression in elevated [CO 2 ] including genes in the TCA cycle and the electron transport chain (Figs 3, 5, S3, S4 and Table S4), also supporting recent findings for dark respiration in elevated [CO 2 ] for young leaves (Markelz et al., 2014). In contrast, spring plants of the same developmental leaf age responded differently to elevated [CO 2 ] and in general responded less. Although acclimation of the photosynthetic machinery occurred, in line with what has been observed in FACE experiments, where photosynthesis is often downregulated following long-term exposure to elevated [CO 2 ] (Ainsworth & Long, 2005), for spring plants, this response to elevated [CO 2 ] was far less apparent. Similarly, the response of transcripts associated with the respiratory pathway to elevated [CO 2 ] also differed between spring and control plants, with increased and decreased transcriptomic responses for spring and control plants, respectively. This reprogramming of the respiratory transcriptome is of interest and is likely to be most important in explaining the enhanced growth and biomass production observed in the spring plants irrespective of CO 2 treatment (Fig. 1). Taken together, this study has provided powerful evidence of the link between gene expression changes between spring and control plants that underpin key processes for adaptation. Photosynthetic acclimation to elevated CO 2 in the control plants as often observed in the literature appears absent in the spring plants, suggesting that transcriptional reprogramming of gene loci in the photosynthetic and respiratory machinery is a key target for future studies of adaptation to high CO 2 . In general, spring plant had higher above-ground biomass than control plants but partitioned less biomass to roots (data not shown) and this suggests that the regulation of carbon allocation and control of growth differed fundamentally in these two plant groups. It provides evidence that the control of carbon acquisition, allocation and growth is likely to represent a set of adaptive traits that we have elucidated at the level of gene expression. These results for acclimation and adaptation are summarized in Fig. 5. Plant acclimation and adaptation to future high CO 2 . A schematic overview of plant acclimation and adaptation to e[CO 2 ], considering gene expression, physiological mechanisms and growth responses following short (acclimation)-and long-term (adaptation). For gene expression categories, all annotated genes were averaged, within particular function categories, as defined in MapMan including Calvin cycle, light reactions of photosynthesis, TCA and mitochondrial electron transport, cell wall and secondary metabolism. Full statistical analysis of pertinent gene categories is given in Tables S3 and S4, for agriGO analysis in Figs S3 and S4. Growth and morphological responses were quantified in the controlled environment experiment described in Materials and Methods, with statistical analysis in Table S3.
Large gene expression changes imply that the two populations are adapted to their sites of origin at the level of gene expression. If the two populations were locally adapted, one might expect that genetic differentiation at the sequence level had accumulated between the populations due to a reduction in gene flow because of immigrant nonviability. The evidence here suggests that there is little genetic differentiation at the mRNA sequence level underlying local adaptation, which stands in stark contrast to the results from the DE analysis. Taking both polymorphism and F ST into account, there are indications that the spring population is derived from the control population, and there is evidence that the spring population has experienced a weak but significant genetic bottleneck. To our knowledge, these are the first population genetic data to be reported for plants following long-term exposure to elevated [CO 2 ].

DE of stomatal patterning genes may underlie adaptive phenotypes in Plantago
The role of stomatal numbers as an adaptive trait is questionable although there is no doubt that stomatal numbers are sensitive to atmospheric [CO 2 ] over geological time and also to other aspects of the environment including drought. Thus, it is imperative to explore this phenomenon in greater detail. Because of the increased stomatal numbers in the spring plants in response to elevated [CO 2 ] and the absence of this difference in the control plants, we reasoned that gene expression differences in the well-characterized stomatal patterning pathway (Bergmann & Sack, 2007), known to be responsive to [CO 2 ] (Gray et al., 2000;Hu et al., 2010), might be contributing to the underlying control of this stomatal phenotype. Here, homologues of 13 loci (of which seven were significantly DE) that regulate the initiation and control of symmetrical and asymmetrical cell divisions that are necessary for the formation of leaf stomatal pores in an ordered and regular manner (Bergmann & Sack, 2007;Hu et al., 2010) were identified in the Plantago transcriptome (Fig. 4). This is the first time that stomatal patterning genes further downstream of BCA1 (Engineer et al., 2014) have been implicated in stomatal response to elevated [CO 2 ] and may be the key points of control for the 25-30% of plant species that show increased stomatal production in future elevated CO 2 conditions, rather than the more common patterns of stomatal number decline (Woodward & Kelly, 1995). In Arabidopsis, SCRM2 (Kanaoka et al., 2008) and CDKB1;1 (Boudolf et al., 2004) are both positive regulators of stomatal patterning: gain-of-function SCRM2 mutants exhibit an over-proliferation of stomata, and CDKB1;1 is necessary for meristemoid division and inhibition of satellite meristemoid formation. Thus, the increased expression in the spring plants of both of these genes (Fig. 4) is predicted to give rise to an increase in SI, such as that observed here. Expression of CDKB1;1 is also negatively correlated with ECS (Boudolf et al., 2004), and we found a similar correlation, greater expression of a CDKB1;1 orthologue in the spring plants grown under elevated [CO 2 ] (Fig. 1), highlighting the importance of this gene for the spring plants' locally adapted stomatal phenotype. The physiological consequences of this genomic-scale change for the control of stomatal patterning could be either increased or decreased leaf water loss, as explained above when considering stomatal clumping and diffusion shells. Reduced leaf water loss would appear advantageous leading to improved water-use efficiency, as is often observed in FACE and other experiments, but this stomatal phenotype is less frequent and thus remains to be investigated. In a meta-analysis of FACE experiments, Ainsworth and Rogers (2007) reported no overall significant effect of elevated CO 2 on stomatal numbers, supporting the findings of Reid et al. (2003). These reports suggest that stomatal patterning in response to future [CO 2 ] varies with species and environment and that our understanding of the underlying mechanisms remains to be elucidated. Here, we provide clear evidence that increased stomatal pore number in response to future [CO 2 ] can occur and that this involves critical steps in the regulatory pathways that are known to control stomatal patterning in the model species Arabidopsis (Vat en and Bergmann, 2012).
Plant acclimation and adaptation to a future high CO 2 world?
This study, using a long-term naturally elevated CO 2 spring, has enabled differentiation between plant acclimatory responses and those potentially due to evolutionary adaptation to high CO 2 atmospheres (Fig. 5). We see the signal of plant adaptation to increasing CO 2 , at the level of plant growth and also gene expression (Fig. 5). Key functional attributes that appear adaptive include increased photosynthetic capacity of spring plants and more mitochondrial respiration (Fig. 3). Differential gene expression was considerable between the spring and control plants, with >800 genes DE. Further investigation of these genes using functional analysis in AGRIGO (Figs S3, S4) and MAPMAN (Fig 5) showed that differential gene expression between spring and control plants included cell wall-associated genes such as XTH21, known to influence cell wall loosening and growth and linked to increased cell production in spring plants. Gene expression linked to secondary plant metabolism pathways including flavonoids, terpenoids and phenypropanoids also differed, categories often involved in response to stress, but also acting as increased carbon pools in spring plants. Several categories related to chlorophyll and photosynthesis were also overrepresented, signifying increased gene expression related to enhanced light harvesting, Calvin cycle and respiration in spring plants. These gene expression differences from transgenerational impacts of exposure to elevated CO 2 are the first of their kind and suggest that the spring plants adapted to high [CO 2 ] had increased photosynthetic carbon fixation and that this carbon was used in respiration and to enhance secondary carbon pools as well as to enhance growth. They provide an intruiging insight into the adaptive phenotype in future CO 2 worthy of further study. An intriguing adaptation in stomatal numbers, with increased numbers in spring plants especially in elevated CO 2 , correlated with regulatory steps in the initiation and development of stomatal pores, where the location of origin was significant in determining the expression of key genes controlling stomatal patterning.
In contrast to these data, which strongly suggest [CO 2 ] acting as a selective agent driving evolution, we observed essentially zero differentiation between these two populations at the sequence level. Our comparisons were necessarily limited to the expressed portion of the genome due to the nature of RNA-Seq, and so it is possible that selection has acted recently on standing genetic variation in a small number of regulatory elements outside of coding regions. However, our evidence for very low differentiation is not likely to be down to the small sample size or because we were only able to investigate the mRNA. Indeed, the nucleotide diversity we resolved was at least as high as that in other wild plant populations (Ramos-Onsins et al., 2004;Ingvarsson, 2005;Liu & Burke, 2006). Changes in the expression of key transcription factors and/or microRNAs can impact multiple genes, as with the microRNAs that regulate temperature and CO 2 responses in Arabidopsis (May et al., 2013). Alternatively, adaptation at the gene expression level could be predominantly epigenetic (e.g. inherited parental effects), as reported recently to explain intraspecific variation in relation to environment (Medrona et al. 2014). Future work should examine the role of the epigenome in response to elevated [CO 2 ].

Conclusions
In conclusion, we present the first study of plant response to multigenerational elevated [CO 2 ] characterized using RNA-Seq. We have revealed important functional genomics responses underpinning both adaptation and acclimation to future [CO 2 ] and have found that adaptation over multiple generations to elevated [CO 2 ] has resulted in changes in gene expression without detectable genetic divergence. This strongly suggests a role for regulatory-level adaptation to future elevated [CO 2 ]. Our study has confirmed that some species develop more stomata in response to elevated [CO 2 ], and we have identified two likely key regulators of stomatal patterning that are central in explaining this phenotype. Figure S3. AgriGo analysis of genes differentially expressed at ambient [CO 2 ]. Figure S4. AgriGo analysis of genes differentially expressed at elevated [CO 2 ]. Table S1. Detailed explanation of plant material used in this study. Table S2. Statistical analysis of morphological and growth data from P. lanceolata .  Table S3. Summary of genes which are displayed in Fig. 4. Table S4. A list of all components assembled from the RNA-Seq data.