Co-ordination in Morphological Leaf Traits of Early Diverging Angiosperms Is Maintained Following Exposure to Experimental Palaeo-atmospheric Conditions of Sub-ambient O2 and Elevated CO2

In order to be successful in a given environment a plant should invest in a vein network and stomatal distribution that ensures balance between both water supply and demand. Vein density (Dv) and stomatal density (SD) have been shown to be strongly positively correlated in response to a range of environmental variables in more recently evolved plant species, but the extent of this relationship has not been confirmed in earlier diverging plant lineages. In order to examine the effect of a changing atmosphere on the relationship between Dv and SD, five early-diverging plant species representing two different reproductive plant grades were grown for 7 months in a palaeo-treatment comprising an O2:CO2 ratio that has occurred multiple times throughout plant evolutionary history. Results show a range of species-specific Dv and SD responses to the palaeo-treatment, however, we show that the strong relationship between Dv and SD under modern ambient atmospheric composition is maintained following exposure to the palaeo-treatment. This suggests strong inter-specific co-ordination between vein and stomatal traits for our study species even under relatively extreme environmental change. This co-ordination supports existing plant function proxies that use the distance between vein endings and stomata (Dm) to infer plant palaeo-physiology.


INTRODUCTION
Global diversity in plant and leaf architecture reflects a plasticity in morphology that allows plants to survive in a range of environments (Díaz et al., 2016). In this current era of rapid climate change, understanding the relationships between plant morphological traits and how they might be influenced by the surrounding environment is of the utmost importance, enabling predictions of plant responses over the coming decades as atmospheric carbon dioxide (CO 2 ) rises. Plants are a critical component of the hydrological cycle, influencing the amount of water vapor that is returned to the atmosphere via the process of transpiration (Rodell et al., 2015). The predicted future increases in CO 2 and global temperatures will have an impact on plant physiological function and morphological traits and will consequently influence the hydrological cycle (Gedney et al., 2006;Betts et al., 2007). The present study focuses on vein and stomatal density (SD), two plant morphological traits that play a pivotal role in the transpirational pathway, and attempts to understand how one may influence the other as a plant encounters environmental change.
Stomata are microscopic pores on a leaf surface that regulate gas exchange. CO 2 from the atmosphere which is essential for photosynthesis is exchanged for water vapor from the inside of the leaf (Jones, 1992). Stomata respond to environmental cues, opening in response to increasing light, low carbon dioxide, and high humidity (Assmann, 1999;Outlaw, 2003; See review by Roelfsema and Hedrich, 2005;Vavasseur and Raghavendra, 2005;Shimazaki et al., 2007;Lawson, 2009). Stomatal opening results in an increase in stomatal pore aperture which leads to an increase in both carbon uptake and water loss from the leaf. SD is the number of stomata per mm 2 of leaf tissue and it is determined by various genetic (Nadeau and Sack, 2002;Shpak et al., 2005) and environmental factors (McElwain and Chaloner, 1995;Woodward and Kelly, 1995;Casson and Gray, 2008). A change in SD alters gas exchange along the plants' diffusional pathway, influencing transpiration and therefore water demand. Veins are found in the leaves of plants, and are differentiations of the vascular bundles that transport water and nutrients from the soil to leaves, as well as sucrose from leaves to the storage sites of the plant (Sack and Holbrook, 2006). A network of major and minor veins (some species only have major veins) carries water throughout the leaf tissue to the stomata where it is lost to the atmosphere as water vapor. Vein density (D v ) is the length of veins per leaf area (mm mm −2 ), and in angiosperms it is determined predominantly by the minor veins, as they make up >80% of the total vein length of the leaf (Sack et al., 2012). Minor vein density has been shown to be an important functional plant trait, exerting a strong influence over xylem conductivity (K x ) and outside xylem conductivity (K ox ), parameters that determine leaf hydraulic conductance (K leaf ) (Sack and Frole, 2006;McKown et al., 2010). Thus it could be said that in the same way that SD and size influence the water demands of a plant, the vein architecture influences its water supply.
In order to be successful in a given environment a plant should invest in a vein network and stomatal distribution that ensures balance between both water supply and demand. Maintaining this balance via co-ordinated shifts in venation and stomatal traits should ensure that the plant is operating at optimal efficiency in terms of carbon uptake and water loss, conforming to the optimality principle (Sack and Scoffoni, 2013). Previous studies have found a strong relationship between D v and SD in response to light (Brodribb and Jordan, 2011) and vapor pressure deficit (Carins Murphy et al., 2014) in certain derived plant species and between SD and transpiration (T) across a range of ferns, conifers, and angiosperms from both tropical and temperate ecosystems (Boyce et al., 2009). D v and SD have been shown to be strongly correlated with modeled maximum theoretical stomatal conductance (g max ) in a diverse range of Proteaceae species . Furthermore, a recent study using a range of modern and basal plant species grown in greenhouse conditions has also reported a strong correlation between D v and g max (McElwain et al., 2016b), suggesting that this balance does indeed exist. Moreover, other studies combine anatomical and physiological measurements to uncover the links between the architectural properties of a leaf and photosynthetic potential. For example, a proxy for leaf photosynthetic capacity has been developed based on the mesophyll path length (D m ) between vein endings and stomata. In multi-veined species, veins should be optimally placed to minimize D m ensuring maximum photosynthetic capacity (Brodribb et al., 2007). These results together demonstrate the potential link between leaf hydraulic morphology and photosynthetic physiology and also highlight the ability of plants to maintain a balance between leaf phenotypic traits under environmental change.
The co-ordination between water supply (D v ) and demand (SD) traits is critical to plant success and the relationship between the two seems to be conserved across the major plant groups under present day atmospheric conditions (Boyce et al., 2009;Brodribb and Jordan, 2011;Brodribb et al., 2013;Carins Murphy et al., 2014;McElwain et al., 2016b). However, it is not known whether this relationship is maintained when levels of oxygen (O 2 ) and CO 2 in the atmosphere change. Co-ordination between these two morphological traits may have been critical throughout the past 400 million years during times of fluctuating atmospheric O 2 and CO 2 . Maintaining this balance between water supply and demand may have allowed certain species to operate more efficiently in their environment. For example, it has been widely proposed that the co-evolution of leaf traits (an increase in D v and SD) during the Cretaceous decline in atmospheric CO 2 allowed angiosperms to outcompete other plant groups as they transitioned from predominantly moist to drier habitats (Boyce and Zwieniecki, 2012;de Boer et al., 2012;McElwain et al., 2016b).
Stomatal density has been shown to be inversely proportional to atmospheric CO 2 (McElwain and Chaloner, 1995;Woodward and Kelly, 1995;Beerling and Woodward, 1996;Royer, 2001;Konrad et al., 2008;Franks and Beerling, 2009a,b) and it has been accepted as a palaeo-environmental proxy for CO 2 on this basis. Studies examining SD responses to concurrent changes in atmospheric O 2 and CO 2 are scarce, as are those investigating D v responses to atmospheric change. In one of the few studies, examining both living and herbarium material of Acer monspessulanum L. and Quercus petraea Liebl, no change was observed in D v in response to an increase in CO 2 from 280 to 350 ppm (Uhl and Mosbrugger, 1999). However, in other studies D v has been shown to respond to environmental change (Brodribb and Jordan, 2011;Sack and Scoffoni, 2013;Carins Murphy et al., 2014) and has also been used in models to predict both atmospheric carbon dioxide partial pressure and temperature (Blonder and Enquist, 2014). Furthermore, D v has the potential to be a useful palaeo-environmental proxy, as venation networks are often preserved in fossilized plant material. For example, studies have shown an increase in D v in angiosperms during the Cretaceous period when CO 2 was declining (Feild et al., 2011a;Boyce and Zwieniecki, 2012).
Using a range of early diverging plant species (Figure 1), this study examines for the first time the effect of changing atmospheric conditions on the relationship between D v , SD, and g max . Using four angiosperm and one fern species, the plasticity and co-ordination of these morphological plant traits was assessed in a low O 2 /high CO 2 atmosphere. In the context of this study, co-ordination in plant traits refers to either inter-specific or intra-specific co-ordination. Inter-specific coordination is taken to mean an observable trend in morphological plant traits across the experimental species. Intra-specific coordination on the other hand, is when the direction of change of plant traits is the same within a single species. We acknowledge that the number of species studied here is relatively small and experimental conditions limited to one palaeo-treatment, and as such, results discussed here are merely intended to be a suggestion of the possible behavior of early diverging plant species under changing atmospheric conditions. Further studies using a wider range of species and palaeo-treatments will build on the current study and will allow more robust conclusions to be drawn.
Results were used to determine the robustness of plant function proxies that rely on co-ordination in morphological traits (such as the use of D m as a proxy for palaeo-assimilation rate (Brodribb et al., 2007;Wilson et al., 2015)), specifically when applied at times of the geological past when the atmospheric composition was different from that of today. Theoretical maximum stomatal conductance (g max ) is a plant functional trait calculated using both SD and anatomical measurements of stomatal geometry (Franks and Beerling, 2009a). This trait has FIGURE 1 | Diagram illustrating the phylogenetic placement of selected experimental plant species. Constructed using http://www.mobot.org/MOBOT/research/APweb/. been used in palaeo-CO 2 proxy models (Franks et al., 2014) to infer past CO 2 levels from the stomatal conductance (g s ) of ancient fossil species. The extent of the relationship between g max and operational stomatal conductance (g op ) has been shown in two angiosperm species (Franks et al., 2009;Dow et al., 2014) and more recently in a range of basal angiosperm, gymnosperm, and fern species (McElwain et al., 2016b). Across this range of basal species g max and g op were found to be strongly related (r 2 = 0.54), with a scaling relationship of g op = 0.25 g max . Therefore g max was used in the current study as a means to infer changes in physiological behavior with a change in the concentration of O 2 and CO 2 and to relate this to changes in D v .

Species and Growth Conditions
Plant species from two different evolutionary plant groups, ferns (Cyathea australis) and angiosperms (Chimonanthus praecox, Magnolia delavayi, Cornus capitata, and Zantedeschia aethiopica), were grown for approximately 7 months in Conviron (Winnipeg, MB, Canada) BDW-40 walk-in plant growth chambers at PÉAC (Programme for Experimental Atmospheres and Climate), Rosemount Environmental Research Station, University College Dublin. For each of these plant groups, the earliest diverging species obtainable from each plant family was used in order to follow the nearest living relative (NLR) protocol (Mosbrugger, 2009), whereby the responses of extant plant species can be said to reflect the responses of their extinct relatives. Growth conditions were set to ambient (two chambers at 21% O 2 and 400 ppm CO 2, O 2 :CO 2 ratio of 525) and a palaeo-treatment of low O 2 /high CO 2 (three chambers at 16% O 2 and 1900 ppm CO 2, O 2 :CO 2 ratio of 84.21). These conditions represent prehistoric modeled atmospheres (Bergman et al., 2004;Berner, 2009) that likely occurred multiple times throughout the last 400 million years, for example in the Devonian (∼ 419-359 mya), the late Triassic (∼ 218-201 mya), and Jurassic periods (∼ 201-145 mya) (Willis and McElwain, 2014). Plants were given water and nutrients according to the individual species requirements (See Supplementary Table S1). Chambers were set to a 16 h day/8 h night schedule, with day/night temperatures of 20 • C/15 • C, relative humidity of 65% and a photosynthetic photon flux density (PPFD) of 600 µmol m −2 s −1 .
SD, D v , and g max Quantification SD and D v were quantified on three leaves per plant, and three to four plants per species per treatment using a modified vein density protocol (Berlyn and Miksche, 1976;Perez-Harguindequy et al., 2013). Leaves were cleared in 5% NaOH, bleached, and then brought through a series of 30, 50, 70, and 100% ethanol. Leaves were then stained using Safranin and Fast Green before being brought back through the ethanol series in the reverse order (100-30% ethanol). Leaves were then suspended in distilled water before being mounted on glass slides for microscopy. For SD quantification, multiple images (on average six images) were taken over an area of approximately 1 cm 2 per leaf and stomatal counts performed using Image J software inside a superimposed grid of 0.09 mm 2 on each image (this area was already determined to be the most representative of the entire leaf using the protocol of Poole and Kürschner (1999) from Jones and Rowe (1999)). For D v quantification, images were taken on three areas of each leaf (an area of approximately 2 mm 2 near the tip, center, and bottom of the leaf near the petiole) and using Image J software the length of veins in these areas was manually traced (using a Wacom Intuos4 pen tablet) and D v calculated in Excel. Minor D v is independent of leaf size and accounts for the majority (>80%) of total vein length per area in most angiosperms (Sack et al., 2012), therefore vein length per area was calculated only on minor veins (quaternary orders upward), major veins being excluded from analysis. Images were taken using a Leica DM2500 microscope with Leica DFC300FX camera (Leica R Microsystems, Wetzlar, Germany) attached and Syncroscopy Automontage (Synchroscopy, Cambridge, UK) digital imaging software was used to impose grids and scale bars on each image.
For g max quantification, anatomical measurements of 90 to 120 stomata per species and per treatment were obtained using the same images used for SD determination (See Table 1 for parameter values). g max was calculated using the following diffusion equation (Parlange and Waggoner, 1970;Franks and Beerling, 2009a): (1) dw = diffusivity of water vapor at 25 • C (0.0000249 m 2 s −1 ), ν = molar volume of air (0.0224 m 3 mol −1 ), SD = Stomatal density (m −2 ), pa max = max stomatal pore area (m 2 ), pd = stomatal pore depth (m). The maximum stomatal pore area was calculated (treating the pore as an ellipse) by using stomatal pore length as the long axis, pore length/2 as the short axis and taking the stomatal pore depth as being equal to the width of a fully turgid guard cell (Franks and Beerling, 2009a,b). It is important to note here that guard cells examined in the current study were not experimentally maintained at maximum turgor before anatomical measurements were made. Even though this might lead to a slight underestimation of g max , this approach is in line with that used in many other palaeo-studies (Franks and Beerling, 2009a,b).

Vein density responses of a selection of non-angiosperm species from 2009 palaeo-experiment
In order to assess whether other non-angiosperm species show a change in D v under different atmospheric compositions, dried plant material from a previous palaeo-experiment was analyzed for D v . Representatives of both the gymnosperms (Agathis australis, Lepidozamia peroffskyana, and Ginkgo biloba) and ferns (Osmunda regalis) were grown for 18 months in walk-in Conviron growth chambers in both an ambient and low O 2 /high CO 2 atmosphere (Ambient treatment: 20.9% O 2 and 380 ppm CO 2, low O 2 /high CO 2 treatment: 13% O 2 and 1500 ppm CO 2,

Statistical Analysis
Data were first checked for normal distribution and Generalized Linear Models were run in Minitab (version 16.1.1) statistical software to investigate differences in SD, D v , and g max between treatments. Minitab (version 16.1.1) statistical software was also used for correlation tests, boxplot representation of data, and to graphically display percent changes in SD, D v , and g max . RStudio (version 0.99.489) was used for Standardised Major Axis (SMA) regression analysis and for scatterplot representation of data.

RESULTS
Species show a varied and species-specific response in SD, D v , and g max to the palaeo-treatment (Figure 2; Table 1). A significant decrease in SD is seen in two species (Chimonanthus praecox shows a 24% and Zantedeschia aethiopica a 34% decrease), Magnolia delavayi and Cyathea australis show a non-significant yet noticeable decrease (8 and 15%, respectively), and another angiosperm (Cornus capitata) shows a small (9%), but nonsignificant increase (Figures 2A and 3). D v shows a similar mix of responses, three of the species show a significant decrease (Chimonanthus praecox, Zantedeschia aethiopica, and Cyathea australis decrease by 15, 12, and 17%, respectively), Magnolia delavayi shows a very slight (2%) yet non-significant increase, and Cornus capitata shows a significant (14%) increase (Figures 2B  and 3). Two of the angiosperm species (Chimonanthus praecox and Zantedeschia aethiopica) show a significant decrease in g max (16 and 32%, respectively) in response to the palaeo-treatment, with Magnolia delavayi showing a non-significant decrease (11%), Cornus capitata a non-significant increase (15%), and Cyathea australis a non-significant decrease (2%) (Figures 2C  and 3). See Supplementary Table S2 for results of generalized linear models (F and associated p-values).
Three out of four angiosperm species (Chimonanthus praecox, Cornus capitata, and Zantedeschia aethiopica) and the fern species Cyathea australis demonstrate intra-specific co-ordination of D v , SD, and g max in response to the palaeo-treatment (Figure 3). This co-ordination is evident even though one angiosperm species (Cornus capitata) shows an increase in all three parameters while the remaining species show a decrease. Magnolia delavayi shows co-ordination in two plant traits in response to the palaeotreatment (SD and g max ), however, co-ordination is lacking between both of these parameters and D v .
The positive relationship between D v and SD ( Figure 4A) is strong under ambient conditions (Pearson's correlation coefficient: r = 0.91, SMA regression: r 2 = 0.82) and it persists in the palaeo-treatment (r = 0.86, r 2 = 0.73). Similarly, D v and g max (Figure 4B) show a strong relationship under both the ambient (r = 0.87, r 2 = 0.76) and palaeo-treatment (r = 0.73, r 2 = 0.53). The slopes of the regression lines between D v and SD and D v and g max are not significantly different in the ambient and palaeo-treatments (D v = 0.01SD + 2.53 for ambient and D v = 0.01SD + 2.39 for palaeo-treatment; D v = 0.003g max + 2.85 for ambient and D v = 0.003g max + 2.66 for palaeo-treatment.

DISCUSSION
SD, D v , and g max Responses to Low O 2 /High CO 2 Results of the current study reflect variability in SD responses to atmospheric change. The inverse relationship between SD and CO 2 has been well documented in the literature using both fossil, herbarium, and living plant material (Woodward and Kelly, 1995;McElwain and Chaloner, 1995;Beerling and Woodward, 1996;Royer, 2001;Konrad et al., 2008;Franks and Beerling, 2009a,b). However, this inverse relationship is not universal across all species (Beerling and Kelly, 1997;Haworth et al., 2013). Stomatal responses to O 2 are not well documented in the literature to date. The few existing studies, however, show an increase in stomatal index (ratio of stomata to epidermal cells or SI) in response to growth in 35% O 2 (Beerling et al., 1998), and a range of SD and SI responses to growth in a combined low O 2 /high CO 2 treatment, as well as to separate low O 2 and high CO 2 treatments (Haworth et al., 2013).
The observed decrease in D v in the majority of species exposed to the palaeo-treatment is likely a consequence of an overall lower water demand due to stomatal optimisation in a high CO 2 atmosphere. Reduced g s in response to high CO 2 has been shown in previous studies (Haworth et al., 2013). This overall reduction in SD and D v reflects a balance between water supply and demand in the palaeo-treatment, and the overall result would likely be a reduction in allocation of resources to nonessential veins and stomata, and a maximization of resource use. It is important to acknowledge that the SD and D v responses observed in the current study cannot be attributed to either low O 2 or high CO 2 alone without undertaking additional and separate sub-ambient O 2 and elevated CO 2 growth experiments using the same species. For the purposes of this analysis, suffice it to say that any SD and D v responses are the result of the specific O 2 :CO 2 ratio in the palaeo-treatment growth chambers.
It is noteworthy that only two of the studied species (Chimonanthus praecox and Magnolia delavayi) have vein densities higher than 6 mm mm −2 (Figure 2B), the 'critical vein density' (Feild et al., 2011a;de Boer et al., 2012) that angiosperms are thought to have surpassed as they rose to dominance in the Cretaceous. This emphasizes the similarity between our chosen study species and those very early evolving angiosperms that had vein densities as low as non-angiosperms (Brodribb and Feild, 2010;Boyce and Zwieniecki, 2012). The significant D v decrease seen in the fern species (Cyathea australis) is interesting (Figure 2B), as it is thought that non-angiosperm species are incapable of altering their vein architecture in the same way that angiosperm species can (Boyce et al., 2009;de Boer et al., 2012). Non-angiosperms seem to exhibit limited plasticity in D v when exposed to a long-term palaeo-atmospheric treatment (Figure 5).
Examination of archived leaf material of three gymnosperms and one fern species from a 2009 palaeo-experiment (Haworth et al., 2013) shows that growth in a low O 2 /high CO 2 atmosphere results in a change in D v in one gymnosperm species (Agathis australis), but has no effect on D v in the remaining species (two gymnosperms and one fern). These non-angiosperms have vein densities below 2 mm mm −2 and all have either parallel or dichotomously branching major vein networks, implying that this vein configuration may lack developmental plasticity. The major veins of non-angiosperms are generally thicker in diameter and their xylem anatomy is distinct from that of the angiosperms, lacking an important feature that is believed to be paramount in the proliferation of minor veins, vessels with simple perforation plates. Only angiosperms evolved these less resistive perforation plates, and this in conjunction with the development of thinner minor veins may have allowed this plant group to outperform non-angiosperms (Feild and Brodribb, 2013). Angiosperms also possess vein endings that are diffuse or dispersed throughout the leaf allowing them to develop more reticulate venation patterns, whereas gymnosperms with their marginal vein endings lack this ability (Boyce, 2005). An important implication of the current findings is that some angiosperm species are able to alter their vein density on a developmental time-scale in response to a change in atmospheric composition; studies to date have only discussed CO 2 -driven D v changes across evolutionary timescales (Brodribb and Feild, 2010;Boyce and Zwieniecki, 2012;McElwain et al., 2016b). Results of the current study indicate that at least in some early diverging species, D v is a plant functional trait that can respond dynamically to atmospheric change.
Relationship between D v , SD, and g max The strong relationship observed in both the ambient and palaeotreatment between D v and SD ( Figure 4A) demonstrates interspecific co-ordination in two morphological plant traits that determine hydraulic supply and demand across different plant lineages and under a changing atmosphere. Furthermore, the robust relationship observed between D v and g max (Figure 4B) demonstrates that morphology has the potential to influence the physiological behavior of these species, via the strong relationship already found between g max and g op (Franks et al., 2009;Dow et al., 2014;McElwain et al., 2016b). Examining the direction of change in the morphological traits for each species it is clear that a high degree of intra-specific co-ordination is also occurring (Figure 3). Three out of four angiosperm species and the fern species show intra-specific co-ordination in D v , SD, and g max in response to the palaeo-treatment. Coordination between traits that determine the water relations (supply and demand) of a plant is critical for its survival. For example, an increase in SD and/or g max would increase the evaporative demands of the plant and without a corresponding increase in D v (to match the increase in water demand with an increase in hydraulic supply) the plant would be maladapted to its environment and would most likely not survive. The opposite scenario would not be as detrimental to plant survival, however, a decrease in SD and/or g max without a corresponding decrease in D v would result in a waste of resources, the construction of veins being costly to the plant (Sack and Scoffoni, 2013). This ability to co-ordinate morphological traits under a changing atmosphere likely occurred throughout plant evolutionary history as the composition of atmospheric O 2 and CO 2 fluctuated, allowing certain plant species to adapt and survive. During the Cretaceous decline in atmospheric CO 2 for example, it is thought that angiosperms were able to increase their gas exchange capacity (thereby increasing photosynthetic rates) by evolving smaller stomata (Franks and Beerling, 2009a), and by increasing both the density of stomata on the leaf surface and the density of veins (Boyce et al., 2009;Brodribb and Feild, 2010;de Boer et al., 2012;McElwain et al., 2016b). Furthermore, angiosperms that surpassed the 'critical vein density' of 6 mm mm −2 were able to out-compete the gymnosperms and ferns in niches with high evapotranspirational demand where an increase in water supply to the leaf would have been necessary for survival (de Boer et al., 2012). Higher vein densities have been suggested to confer a higher capacity for CO 2 uptake and an increased range of g op ; this would explain the ability of angiosperms to expand to such diverse habitats and to outcompete species that are more constrained in their venation and hence gas exchange capacity (McElwain et al., 2016b). A recent study suggests that angiosperms are indeed hydraulically optimized for a diverse range of environments, achieving this by maintaining an equal vein to vein and vein to evaporative surface distance in the leaf (Zwieniecki and Boyce, 2014). Ferns are under-invested hydraulically due to their thin leaves and large vein to vein distances, and while some gymnosperms do approach optimal investment by producing thicker leaves in more water-demanding environments, they are as a group sub-optimal in terms of vein placement (Zwieniecki and Boyce, 2014).
The current study supports these theories by demonstrating a higher degree of plasticity in D v in some early diverging angiosperms in response to a changing O 2 :CO 2 ratio, compared to the studied gymnosperms and ferns (Figure 5). It is interesting, however, that examined non-angiosperm species from the 2009 palaeo-experiment show limited plasticity in D v as well as in SD (see SD results for these non-angiosperm species reported in Haworth et al., 2013). This demonstrates that while these species do not show a high degree of morphological plasticity in response to a changing atmosphere on experimental time-scales comparable to the studied angiosperms, FIGURE 5 | Comparison of the D v responses of gymnosperms, ferns, and angiosperms to growth in low O 2 /high CO 2 conditions. Different letters signify a significant difference between treatments. A, ambient treatment; P, palaeo-treatment. * symbols denote outliers. they do demonstrate similar co-ordination in leaf morphological traits.
Furthermore, these results suggest that angiosperms are not only capable of showing morphological plasticity in response to rising O 2 and declining CO 2 (Boyce et al., 2009;Feild et al., 2011a;de Boer et al., 2012), but also to declining O 2 and high CO 2 conditions. The robust positive relationship observed here between the density of veins and stomata strongly supports the theory suggested by Brodribb et al. (2007), whereby multiveined leaves optimize the placement of veins in relation to stomata so that the distance water needs to travel through the resistive mesophyll (D m ) is minimized. Although D m was not directly measured in this study, the morphological coordination observed suggests that any change in D v will elicit a corresponding change in SD or vice versa, allowing the leaf to minimize the distance between veins and stomata, and to maximize photosynthetic performance and operational efficiency of the leaf.

Implications for Past Plant-Atmosphere Interactions
The experimental species examined show a species-specific and varied response to growth in the palaeo-treatment, yet the strong positive relationship between D v and SD persists ( Figure 4A). Furthermore, the positive relationship observed between D v and g max (Figure 4B) demonstrates the link between hydraulic and gas exchange/diffusional processes in these species, as shown in previous studies (Sack et al., 2003;Boyce et al., 2009;Brodribb and Feild, 2010;Feild et al., 2011b;McElwain et al., 2016b). The finding that D m is most likely maintained under changing atmospheric conditions (due to the intra-specific coordination between D v and SD) has important implications when attempting to understand plant-atmosphere interactions throughout the last 400 million years of plant evolution.
A lack of co-ordination in D v and SD on developmental timescales would result in a plant that is morphologically and physiologically out of sync, negatively impacting operational efficiency and overall fitness under changing atmospheric conditions. Furthermore, plant species that exhibited plasticity in these morphological traits under a changing atmosphere would likely have had an ecological advantage over plant species that were morphologically inflexible, being able to maximize their photosynthetic capacity as the surrounding environment changed (McElwain et al., 2016b). The finding that a proxy for photosynthetic capacity (D m ) (Brodribb et al., 2007) remains stable under changing atmospheric conditions is important for accurate initial parameterisation of mechanistically based models used to predict palaeo-CO 2 (Franks et al., 2014) since these require robust estimates of palaeo-assimilation (Franks et al., 2014;McElwain et al., 2016a). The current study focuses on plant responses to an experimentally imposed low O 2 /high CO 2 atmosphere. This O 2 :CO 2 ratio occurred multiple times throughout plant evolutionary history based on model and proxy estimates (Royer, 2001;Bergman et al., 2004;Royer et al., 2004;McElwain et al., 2005;Berner, 2006Berner, , 2009Steinthorsdottir et al., 2016), however, it would be beneficial to investigate the effect of other atmospheric O 2 :CO 2 ratios on these plant trait relationships in order to test their linearity. Further experiments examining plant responses to a range of palaeo-atmospheric conditions will build on these results, providing a better picture of plant-atmosphere interactions over the past 400 million years and allowing predictions of future plant responses to global climate change.

CONCLUSION
Species show a varied response in SD, D v , and g max to growth in an experimental low O 2 /high CO 2 palaeo-atmosphere.
Regardless of this variation in responses, a strong relationship is observed between D v and SD and D v and g max under both the ambient and palaeo-atmosphere. Gymnosperms studied here appear to lack the same degree of developmental plasticity in D v compared to the angiosperms, at least on short experimental time-scales. The ability to increase their range of D v values may have contributed to the success of angiosperms during the Cretaceous decline in CO 2 ; a high degree of plasticity in this trait possibly provided early diverging angiosperms with a competitive advantage over other seed plant groups in more changeable environments. The tight relationship observed between D v and SD in the palaeo-treatment suggests that D m is likely maintained under environmental change and lends confidence to existing palaeo-CO 2 proxies that use this parameter in their models. Further studies examining the robustness of these plant trait relationships under a range of O 2 :CO 2 ratios are needed in order to elucidate the full spectrum of plant-atmosphere interactions throughout the last 400 million years.

AUTHOR CONTRIBUTIONS
CE-F carried out vein density and SD analysis, anatomical stomatal measurements for g max calculation, statistical analysis, and drafted the manuscript. AP contributed SD and g max data for Cyathea australis. CE-K was involved in the 2009 palaeoexperiment and therefore provided dried leaf material for vein density analysis. JM is the principal investigator. All authors read, revised, and approved the final manuscript.

FUNDING
We gratefully acknowledge funding from a European Research Council grant (ERC-279962-OXYEVOL).

ACKNOWLEDGMENTS
We thank Mr. G. Kavanagh and Ms. B. Moran (UCD, Ireland) for their technical assistance. We would also like to thank both J. D. Fitz. Gerald and Dr. S. P. Batke for their assistance with statistical analysis. Finally, we thank the Editor and both reviewers for their valuable comments and suggestions.