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₂ and 400 ppm CO_2, O₂:CO₂ ratio of 525) and a palaeo-treatment of low O₂/high CO₂ (three chambers at 16% O₂ and 1900 ppm CO_2, O₂:CO₂ 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 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² per leaf and stomatal counts performed using Image J software inside a superimposed grid of 0.09 mm² 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² 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^® 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):

dw = diffusivity of water vapor at 25°C (0.0000249 m² s^-1), ν = molar volume of air (0.0224 m³ mol^-1), SD = Stomatal density (m^-2), pa_max = max stomatal pore area (m²), 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).

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₂/high CO₂ atmosphere (Ambient treatment: 20.9% O₂ and 380 ppm CO_2, low O₂/high CO₂ treatment: 13% O₂ and 1500 ppm CO_2, O₂:CO₂ ratios 552.63 and 86.67, respectively). These species have either parallel or dichotomously branching major veins, therefore it was not necessary to clear and stain leaves for D_v observation. D_v was calculated on one leaf per plant and three plants per species per treatment. Leaves were imaged using a Nikon SMZ1000 stereomicroscope with Leica DFC490 camera attached (Leica^® Microsystems, Wetzlar, Germany), and using Image J software the length of veins was manually traced and D_v calculated in Excel in an area between 20 mm² and 60 mm².

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 of the current study reflect variability in SD responses to atmospheric change. The inverse relationship between SD and CO₂ 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₂ 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₂ (Beerling et al., 1998), and a range of SD and SI responses to growth in a combined low O₂/high CO₂ treatment, as well as to separate low O₂ and high CO₂ 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₂ atmosphere. Reduced g_s in response to high CO₂ 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 non-essential 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₂ or high CO₂ alone without undertaking additional and separate sub-ambient O₂ and elevated CO₂ 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₂:CO₂ 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₂/high CO₂ 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₂-driven D_v changes across evolutionary time-scales (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.

The strong relationship observed in both the ambient and palaeo-treatment between D_v and SD (Figure 4A) demonstrates inter-specific 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. Co-ordination 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 mal-adapted 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₂ and CO₂ fluctuated, allowing certain plant species to adapt and survive. During the Cretaceous decline in atmospheric CO₂ 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₂ 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₂:CO₂ 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, 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₂ and declining CO₂ (Boyce et al., 2009; Feild et al., 2011a; de Boer et al., 2012), but also to declining O₂ and high CO₂ 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 multi-veined 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 co-ordination 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.

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 co-ordination 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 time-scales 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₂ (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₂/high CO₂ atmosphere. This O₂:CO₂ 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, 2006, 2009; Steinthorsdottir et al., 2016), however, it would be beneficial to investigate the effect of other atmospheric O₂:CO₂ 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.