Extremely thick cell walls and low mesophyll conductance: welcome to the world of ancient living!

Highlight Low mesophyll conductance was the most important limiter of photosynthesis in the majority of evolutionarily old species and resulted from thick cell walls and low exposed chloroplast area.


Introduction
The 36-fold variation of net photosynthetic rate across C3 species (Wright et al., 2004) cannot be explained only by stomatal restrictions and biochemical potentials, because the CO 2 diffusion efficiency from substomatal cavities to chloroplasts (mesophyll diffusion conductance; g m ) plays an important role in shaping photosynthetic capacities and leaf resourceuse efficiency across Earth's ecosystems (Warren et al., 2003;Niinemets et al., 2009bNiinemets et al., , 2011Terashima et al., 2011;Tosens et al., 2012bTosens et al., , 2016. The general consensus is that ignoring g m in carbon gain calculations prevents the correct simulation of photosynthesis due to biased maximum velocity of carboxylation especially in stressful field environments (Manter and Kerrigan, 2004;Niinemets et al., 2011). Still, g m has been ignored in carbon gain calculations because of our lack of knowledge about its various causes and the main drivers of its variation across species (Niinemets et al., 2009a).
Mesophyll conductance has now been estimated for more than 100 species from all major plant groups. Available data show that considerable variation of g m between plant groups exists and it is larger than that of stomatal conductance (g s ) Tosens et al., 2016). g m depends on the length and characteristics of the CO 2 diffusion path from substomatal cavities to chloroplasts Hassiotou et al., 2009;Terashima et al., 2011;Tosens et al., 2012b). The main limitations to CO 2 diffusion are in the mesophyll liquid phase. Two traits, mesophyll cell wall thickness (T cwm ) and chloroplast surface area exposed to intercellular airspace (S c /S), have been highlighted as the strongest limiting factors of g m , but these anatomical traits are highly variable among species Terashima et al., 2011;Tosens et al., 2012b). In addition, recent studies suggest that chloroplast thickness is also an important barrier limiting diffusion of CO 2 to Rubisco (Tosens et al., 2012a(Tosens et al., ,b, 2016Peguero-Pina et al., 2012;Tomás et al., 2013). The reduction of photosynthesis due to g m can be about 25% for mesophytic species and up to 75% in evergreen species . The latter have more robust leaves, e.g. thicker cell walls due to their adaptations to their specific environmental stressors (Niinemets, 2016). Contrarily, ferns and allies have seemingly soft leaves and low leaf mass per area (LMA), but the extent of the photosynthetic limitation by g m and net assimilation are in a similar range to that reported for sclerophyllous angiosperms due to their very thick mesophyll cell walls and very low S c /S (Niinemets et al., 2009b,c;Tosens et al., 2016). That is, g m and its underlying traits cannot be simply extended across plant groups, and additional factors such as evolutionary age can play an important role in its determination.
Across spermatophytes g m is lowest in gymnosperms . Six-fold variation in g m has been shown among the 13 conifer species studied so far (Warren et al., 2003;De Lucia et al., 2003;Manter and Kerrigan, 2004;Black et al., 2005;Mullin et al., 2009;Peguero-Pina et al., 2012, 2016. Collectively, the available data suggest that g m is an important limiting factor of photosynthesis across and within conifer species. Despite its undoubted importance, only Peguero-Pina et al. (2012, 2016 have investigated g m together with the underlying structural correlates. Therefore, integrative approaches are needed to characterize the nature of g m and its importance in controlling realized photosynthesis rates in gymnosperms.
The leaf economics spectrum (LES) is a dimension of ecological variation reflecting differences across species in the cost of the investment in a unit of leaf area and the return on that investment (Wright et al., 2004). Despite the broad use of LES relationships in generalizing global leaf structure-function relationships, the morpho-physiological basis underlying given trade-offs are poorly understood, especially for species placed at the low-return end (Shipley et al., 2006;Niinemets et al., 2009a,b). These species are characterized by high LMA and robust structure, and therefore the realized photosynthetic rates are modified in an important way by leaf morphology (Flexas et al., 2008;Hassiotou et al., 2009;Niinemets et al., 2011;Tosens et al., 2012b). For example, high LMA could be associated with high net assimilation rate (A area ) if its main driver is great thickness of mesophyll and high S c /S, or with low A area if it is connected with its underlying anatomical variations such as thick mesophyll cell walls (Hassiotou et al., 2009;Tomás et al., 2013). However, the integrative trait of LMA may fail for species like Australian Proteaceae, which have mesophytic mesophyll tissue embedded in highly sclerophytic tissue, or in ferns and allies, which have low LMA but high T cwm resulting in low photosynthesis (Niinemets et al., 2009c;Tosens et al., 2012bTosens et al., , 2016. Understanding the underlying ultrastructural traits controlling the global LMA-photosynthesis-nitrogen relationships is inescapable for understanding photosynthetic patterns across Earth's ecosystems. The reasons for lower photosynthesis rates and intrinsic water use efficiency (WUE i ) in gymnosperms are currently poorly known due to lack of studies separating mesophyll and stomatal diffusional and mesophyll biochemical photosynthetic limitations. Information about g m with its underlying structural traits is especially limited. Here we set out to characterize the spectrum of photosynthetic strategies and related physiological and morphological traits in gymnosperms positioned at the lower return end of the leaf economics spectrum based on 11 gymnosperm and two evolutionarily older species, Psilotum nudum and Selaginella uncinata, in order to enhance the evolutionary context. Overall, the included species are evolutionarily old with a phylogenetic age extending from 306 million years for Psilotaceae to 75 million years for the genus Podocarpus (Pryer et al., 2004;Biffin et al., 2012). Metasequoia glyptostroboides, Cycas revoluta, Macrozamia riedlei and Araucaria heterophylla can even be regarded as living fossils as their morphology and habitat have a high resemblance to fossils dated to the Mesozoic (>110 My ago) and even to the late Paleozoic (300 My ago) (e.g. Norstog and Nicholls, 1997;Rydin et al., 2004;Zhang et al., 2015). The specific aims of the study were (i) to analyse the photosynthetic limitations to understand if diffusional limitations, especially g m , are the most substantial constraints in evolutionarily old plants; (ii) to assess which structural traits are responsible for low CO 2 diffusion conductance in evolutionarily old species; and (iii) to understand the morpho-physiological basis underlying LES traits in evolutionarily old species.

Plant material and growth conditions
Thirteen evolutionarily old species (Fig. 1 and Supplementary Table  S1 at JXB online) with widely varying shape, size, longevity and structure of photosynthetic organs covering a broad range of evolutionary ages were included in the analysis. Out of the selected species, 11 (Araucaria heterophylla, Cupressus sempervirens, Cycas revoluta, Ephedra minuta, Macrozamia riedlei, Metasequoia glyptostroboides, Picea abies, Pinus sylvestris, Podocarpus alpinus, Podocarpus nivalis and Taxus baccata) were gymnosperms from three of the four gymnosperm divisions. Two primitive herbaceous plant species (a whisk fern, Psilotum nudum, and a clubmoss, Selaginella uncinata) were also included.
Five-month-old saplings of M. glyptostroboides, P. alpinus, P. nivalis, S. uncinata and T. baccata were obtained in October 2014 from a local nursery and transferred to a controlled-conditions phytotron, where they grew for 2-4 months before the measurements were begun. Other species were grown in the phytotron from seeds or seedlings brought from their country of origin (see Supplementary  Table S1). All measured leaves emerged and matured in the phytotron conditions. Relative humidity was maintained at 60% and air temperature at 25 °C/22 °C (day/night). Light was provided by high pressure sodium lamps (400 W, E40 tubular clear, Beijing Luxram Lighting Ltd, China) for a 13 h photoperiod. The photosynthetic photon flux density (Q) incident to the plants was maintained at 1000 µmol m -2 s -1 reflecting the natural growth light conditions of the species (Supplementary Table S1). Plants were watered every second day with tap water and fertilized with soluble fertilizer once a week (Substral Miracle for evergreens, N:P 2 O 5 :K 2 O:MgO 23:9:12:2 with microelements). In all cases, fully mature non-senescent foliage was used for the measurements. Low (P. nudum, left) and high (P. sylvestris, right) chloroplast area exposed to intercellular airspace (S c /S). Arrows in (C) indicate sparsely (left) and tightly packed chloroplasts (right). The P. sylvestris cross-section exhibits lobed cells (circled in C).

Foliage gas exchange measurements and estimation of mesophyll conductance
Simultaneous gas exchange and fluorescence measurements were conducted with a GFS-3000 portable photosynthesis system (Walz GmbH, Germany). Foliage elements were arranged side by side in the leaf chamber to avoid any overlap and a digital photograph was taken for estimation of the exact projected foliage area enclosed in the chamber after its closure. In all measurements, vapor pressure deficit was kept at ca 1.5 kPa, and leaf temperature, measured with a thermocouple, at 25 ºC. Leaf net photosynthesis, the steady-state net photosynthesis rate under saturating Q of 1000 µmol m -2 s -1 (90% red and 10% blue light) and 400 µmol CO 2 mol -1 (atmospheric CO 2 concentration, C a ), was recorded after the stomata opened and leaf gas exchange rates reached maximum steady state, 30-60 min after the enclosure of foliage in the chamber. Net assimilation vs. ambient CO 2 response curves were measured at a Q of 1000 µmol m -2 s -1 by varying C a between 50 and 2000 µmol CO 2 mol -1 . After each steady-state net assimilation rate had been estimated, the steady-state fluorescence level (F) was recorded and a saturating pulse was given by the leaf chamber fluorimeter of the GFS-3000 system (PAM-fluorometer 3055-FL) to measure the maximum fluorescence yield (F m ′) and estimate the effective photosystem II quantum yield (Φ PSII ) as (F m ′-F)/F m ′ (Genty et al., 1989). The area of the leaf chamber was 8 cm 2 . Once the measurements for a CO 2 response curve were completed, the light was switched off and the mitochondrial respiration rate (R n ) was measured after a minimum of 30 min of dark adaptation. g s and intercellular CO 2 concentration (C i ) were calculated according to von Caemmerer and Farquhar (1981).
Mesophyll conductance was calculated according to Harley et al. (1992): where A area is the net assimilation rate, J ETR is the rate of photosynthetic electron transport derived from chlorophyll fluorescence measurements, C i is the CO 2 concentration in sub-stomatal cavities, R d is the rate of non-photorespiratory respiration in the light and Γ* is the hypothetical CO 2 compensation point without R d . R d was estimated as half of R n (Niinemets et al., 2005). This common adaptation is supported by previous experimental observations (Villar et al., 1995;Piel et al., 2002). Γ* was taken as 42.9 μmol mol -1 at 25 °C (Bernacchi et al., 2001) as it has been shown before on vastly variable plant functional groups and g m values that there is no significant difference between g m calculated with Γ* according to Galmés et al. (2005) and Bernacchi et al. (2001) (Tomás et al., 2013). J ETR was estimated from non-photorespiratory conditions (2% oxygen) from an A-C i curve, although the calculations of g m were relatively insensitive to moderate variations in J ETR (<0.5%; see Niinemets et al., 2005). Therefore, the main assumption of Harley's variable J method was not breached. The data were corrected for chamber leaks according to Flexas et al. (2007a) and Rodeghiero et al. (2007) using empty chamber measurements. Leaf absorption was measured with an integrating sphere (leaf absorption coefficient varied between 0.84-0.90). Gaskets were changed regularly, and modelling paste and paraffin film were used with thick leaves to ensure no leakage. Mesophyll conductance was calculated from measurements of net assimilation rate over the C i range of 150-350 μmol mol -1 , because the g m values are stable over this range and its estimates are relatively insensitive to minimal Γ*, R d and A errors (Harley et al., 1992;Niinemets et al., 2006). WUE i was defined as A area /g s . As some plants studied here have small leaves or needles growing on several sides of the branch, not all of these were perpendicular to the light source. However, this is also the case for when they are growing in natural conditions as the branches were positioned in the same direction to the light source as they had grown.
Additionally, the one-dimensional within-leaf gas diffusion model of Niinemets and Reichstein (2003) was applied to provide an anatomically based estimate of g m as in Tosens et al. (2012b) and from A-C i curves as in Ethier and Livingston (2004).

Quantitative limitation analyses of A area and partial limitations of g m
The relative controls on A area imposed by stomatal conductance (l s ), mesophyll conductance (l m ), and biochemical capacity (l b ) were calculated following the Grassi and Magnani (2005) approach, as has been successfully used for a large number of highly variable species from lycophytes, horsetails and ferns to Australian sclerophylls in Tomás et al. (2013) and Tosens et al. (2012bTosens et al. ( , 2016. These three values sum to 100% and characterize the extent to which any of the three limitations curbs photosynthesis at the given values of the other two. The contribution of different components of cellular resistance to total cellular resistance to CO 2 diffusion was estimated from the anatomical model following Tosens et al. (2016). This share of limitation (l i ) by different liquid phase components was calculated as: where l i is the limitation by the cell wall, the plasmalemma, cytosol, chloroplast envelope and stroma, and g i refers to the diffusion conductance of each corresponding diffusion pathway. The limitation of each cellular component was scaled up with S c /S.

Estimation of leaf dry mass per area and nitrogen content
Foliage used for gas exchange measurements was further used for anatomical, morphological and chemical measurements. Individual foliage elements (leaves, needles or scales) were removed from the twigs and scanned at 300 dpi. The foliage was thereafter oven-dried at 70 °C for 48 h, and its dry mass was determined. Leaf projected area was determined from scanned images with ImageJ 1.48v software (Wayne Rasband/NIH, Bethesda, MD, USA) or calculated from the measured length and widths of sides measured with digital precision calipers (Mitutoyo CD-15DC, Mitutoyo Ltd, Andover, UK) for needles and cladodes . From these measurements, LMA and density (D leaf =LMA/T leaf ) were calculated (Niinemets, 1999;Poorter et al., 2009). Nitrogen content per dry mass (N mass ) was determined by the dry combustion method with a Vario MAX CNS elemental analyser (Elementar, Hanau, Germany), and leaf nitrogen content per area as N mass /LMA.

Light and electron microscopy
A sub-sample of leaves used for gas exchange measurements was taken for leaf anatomical measurements. Light and electron microscopy sample preparation and analyses followed Tosens et al. (2012b). Six samples were taken from three specimens of each species. Foliage cuts of approximately 6 × 4 mm (or width of the needle or cladode according to species) were taken from intercostal areas. In C. sempervirens scales and in S. uncinata whole leaves were stripped from the branches with tweezers. The plant material was infiltrated in a fixation buffer (3% glutaric aldehyde and 2% paraformaldehyde in 0.1 M phosphate buffer, pH 6.9) under vacuum in a syringe. The samples were post-fixed for 1 h in an osmium tetroxide solution (2%) and dehydrated in a series of increasingly stronger ethanol solutions and embedded in LR white resin (Electron Microscopy Sciences, Hatfield, PA, USA) according to standard procedures (Tosens et al., 2012b). Subsequently, the samples were polymerized in an oven at 60 °C for 26 h. Semi-thin cross-sections of 1 µm for light microscopy, and ultra-thin crosssections of 70 nm were prepared by an ultramicrotome (Leica EM UC7, Leica Vienna, Austria) for transmission electron microscopy (TEM). The semi-thin sections were stained with toluidine blue for light microscopy, and the sections for TEM were stained with lead citrate and then mounted on formvar-carbon covered copper meshes (Electron Microscopy Sciences). Stained semi-thin sections were viewed in bright field ensuring that all sections were ideally straight with a Nikon Eclipse E600 microscope with phase contrast at magnifications of ×100, ×200 and ×400 and photographed with a Nikon 5 MP digital microscope camera DS-Fi1 (Nikon Corp., Kyoto, Japan). A Philips Tecnai 10 TEM (FEI, Eindhoven, Netherlands) was used to view the ultra-thin sections with the accelerating voltage of 80 kV and magnification between ×1800 and ×14 000. Between five and seven fields of view per sample were measured. From digital images the following parameters were measured: fraction of intercellular airspace (F ias ), leaf thickness from adaxial to abaxial cuticle (T leaf ), mesophyll thickness excluding vascular bundles (T mes ), cytoplasm thickness (T cyt ), chloroplast thickness (T chl ), mesophyll cell wall thickness (T cwm ), chloroplast and mesophyll cell wall area from which S c /S and mesophyll area exposed to intercellular airspace (S m /S) were calculated according to Evans et al. (1994) (Table 2). These characteristics were measured for at least ten spongy and ten palisade parenchyma cells and tissue-volume weighted averages were calculated in the species with distinctive separation of mesophyll to spongy and palisade parenchyma. Thirty cells were analysed for every specimen. Light and electron micrographs were analysed with ImageJ 1.48v software.

Statistical analyses
For each species, all measurements were replicated at least with three individual plants. Linear and non-linear regression analyses, t-tests and ANCOVA were used to examine the relationships among the traits and test for trait differences among species using Statistica 10.0 (StatSoft Inc., Tulsa, OK, USA). The choice between linear and non-linear models depended on the shape of the relationships, degree of explained variance and normality of data and residuals, and models providing the greatest r 2 and lowest deviations from normality of residuals were used. Due to high variability of foliage shapes and structures, individual species were occasionally outliers in statistical relationships, especially in relationships exploring the effects of individual underlying traits on composite variables. For instance, some correlations were significant when the gymnosperms were considered, but not when the spore-bearing species were included. We denote these and analogous outliers separately in bivariate statistical relationships.

Phylogenetic analyses
Sequences of the plant barcode genes maturase K (matK; a gene that encodes an intron splicing protein) and RbcL (a gene that encodes the large subunit of Rubisco; both genes are in the chloroplastic DNA; Hollingsworth et al., 2009) were extracted from NCBI GenBank (www.ncbi.nlm.nih.gov, last accessed July 2015). This information was not available for Macrozamia riedlei, so the closest available relative, M. moorei, was used. Phylogenetic analyses were conducted and phylogenetic trees were created by the neighbor-joining method using standard procedures in MEGA6 (Tamura et al., 2013). The maximum composite likelihood model was used for estimates of evolutionary divergence (Tamura et al., 2004), and the maximum likelihood method based on the Tamura-Nei model (Tamura and Nei, 1993) was used to create phylogenetic topology with the superior log-likelihood value. The evolutionary age of genera was based on previously published literature (Table 1: smaller numbers indicate evolutionarily closer species). This information was correlated with the study parameters using t-tests and ANCOVA. The Akaike information criterion was used as a measure of the relative quality of the multiple linear models used here.

Variation in leaf anatomy and morphology in evolutionarily old species
In order to understand the mesophyll diffusional limitations to photosynthesis, a complete structural and ultrastructural analysis was performed of the photosynthetic organs of all the species. A large variability was observed in their macroscopic anatomy (Fig. 1A), but also structural and ultrastructural parameters (Figs 1B-D and 2), e.g. a P. sylvestris cross-section exhibited lobed cells, increasing mesophyll surface area exposed to the intercellular airspace (Fig. 1C). LMA varied about 7-fold between species, from 45 ± 12 g m -2 in S. uncinata to 308 ± 22 g m -2 in C. sempervirens. Leaf density (D leaf ) and thickness varied 16-and 24-fold, respectively. T cwm varied 5.5-fold (from 0.2 ± 0.05 µm in S. uncinata to 1.2 ± 0.37 µm in P. nivalis), S c /S varied 3.5-fold (from 4.9 ± 0.47 m 2 m -2 in S. uncinata to 17.1 ± 0.62 m 2 m -2 in P. sylvestris), and S m /S varied 2.9-fold (6.1 ± 0.18 to 17.5 ± 0.72 m 2 m -2 ). The variation in LMA was attributed to variation in both of its components, but there was a stronger positive relationship with D leaf (excluding S. uncinata due to its extremely thin leaves) than T leaf (Fig. 3A, B). However, there was no correlation between LMA and mesophyll cell wall thickness (Fig. 3C). There was a strong positive relationship between S c /S and S m /S (see Supplementary Fig. S1A), but neither of them was related to LMA (Fig. 3D for S c /S; r 2 =0.077, P=0.36 for the correlation with S m /S). Similarly, S m /S was not correlated to fraction of intercellular airspaces, F ias (r 2 =0.007; P=0.78). S c /S was unrelated to T mes ( Supplementary Fig. S1B), but S m /S had a positive correlation (Supplementary Fig. S1C). However, this correlation was not significant when only gymnosperms were included (r 2 =0.068; P=0.37).

Interspecific variation in photosynthetic capacity and its relationships with diffusional and structural determinants
Area-and mass-based net assimilation (A mass ) varied about 5-and 8-fold, respectively (from 2.1 ± 0.34 µmol m -2 s -1 in P. nudum to 11.0 ± 0.78 µmol m -2 s -1 in P. sylvestris, and from 21 ± 3.9 nmol m -2 s -1 in C. revoluta to 165 ± 9 nmol m -2 s -1 in M. glyptostroboides). Mesophyll conductance per unit leaf area (g m/area ) varied over 12-fold (from 10 ± 2.1 mmol m -2 s -1 in Ephedra minuta to 124 ± 9 mmol m -2 s -1 in P. sylvestris) and mass-based g m (g m/mass ) over 15-fold (from 0.08 ± 0.005 mmol g -1 s -1 in P. nivalis to 1.20 ± 0.06 mmol g -1 s -1 in M. glyptostroboides). A area depended strongly on g s (Fig. 4A). Importantly, A area scaled positively with g m/area (Fig. 4B), but there was stronger positive correlation between A mass and g m/mass (Fig. 4C), which was significant even with the exclusion of M. glyptostroboides (r 2 =0.42; P=0.016). A similar strength negative curvilinear relationship was found between C a -C i vs g s (Fig. 4D) and the CO 2 drawdown from the intercellular airspace to chloroplasts (C i -C c ) vs g m/area (Fig. 4E). However, a stronger relationship was found between C i -C c vs g m/mass when S. uncinata and M. glyptostroboides were excluded (Fig. 4F).
In addition to large variations in LMA and its components, nitrogen content per dry mass (N mass ) varied 6-fold (from 0.45 ± 0.02% to 2.72 ± 0.07%) and area-based nitrogen (N area ) varied 12-fold (0.6 ± 0.01 to 6.1 ± 0.09 g m -2 ). A mass scaled positively with N mass , but the only deciduous conifer in this study, M. glyptostroboides, diverged from this relationship by having higher A mass (Fig. 5A). There was no relationship between N area and A area (r 2 =0.20; P=0.12). Similarly, A mass was negatively correlated with LMA, while there was no relationship between A area and LMA (Fig. 5B, C).
Analysis of the ultra-anatomical controls on g m demonstrated that g m/area was positively correlated with S c /S. On the other hand, g m/mass and g m/area were negatively associated with T cwm . C i -C c scaled positively with T cwm across gymnosperms, but not across the whole sample (Fig. 6). T cwm did not correlate with density (r 2 =0.030; P=0.55).
In addition to gas exchange-fluorescence methodology, g m was estimated by two alternative methods: it was modeled from mesophyll anatomical traits by using the anatomical model of Niinemets and Reichstein (2003) and from A-C i curves (Ethier and Livingston, 2004, Fig. 7). g m calculated from anatomical measurements correlated better with the estimates obtained from gas exchange-fluorescence measurements than g m obtained from A-C i curves (average discrepancy between g m from gas exchange-fluorescence was 27% with g m from anatomy and 33% with g m from Ethier and Livingston (2004)). Based on the quantitative limitation analysis, the photosynthetic capacity was strongly limited by both g s (range 15-70%) and g m/area (range 12-74%), while limited biochemical capacity (range 8-33%) restricted A area less than stomata and mesophyll (Fig. 8A). With respect to the structural limitations of g m , the estimated gas phase limitation inside the leaf was <1% of the total limitations for all the species (data not shown). Among all the components of liquid phase limitations, cytoplasm and membrane limitations played a minor role, whereas the predominant limitation was exerted by cell walls and chloroplast stroma ( Fig. 8B and Table 2).
The intrinsic water use efficiency (WUE i ) varied ca 3-fold across species. It was negatively correlated with the share of photosynthesis limited by g m/area (Fig. 9A). WUE i also scaled positively with the photosynthesis limitation by stomata (Fig. 9B). When considered together, the interspecific variation in WUE i was driven by both l s and l m (see Supplementary  Table S2).

The effect of divergence time on structure and physiology
In the given dataset, LMA and A area were negatively correlated with the estimated evolutionary age of the genus (Supplementary Fig. S2 and Supplementary Table S3a; r 2 =0.45; p=0.012). However, g m/area itself was not related to plant evolutionary age (r 2 =0.055; p=0.44) and g m/area influenced A area regardless of evolutionary age (Supplementary  Table S3b). Nevertheless, the lowest mesophyll conductances were observed for the oldest genera (Ephedra: 10.3 mmol m -2 s -1 ; Psilotum: 12.5 mmol m -2 s -1 ; Selaginella: 20.5 mmol m -2 s -1 ). However, the highest average g m/area was recorded for Pinus (124 ± 9 mmol m -2 s -1 ), although this genus is older than four others (Picea, Podocarpus, Taxus, Cupressus) in this study. In addition, g s correlated negatively with the age of the genus (r 2 =0.41; P=0.017), but there was no significant correlation of WUE i with the evolutionary age of the genus (r 2 =0.083; P=0.34).

Structural and ultrastructural traits determining LES relationships in evolutionarily old species
LMA in our study depended on leaf density while its relationship with T leaf was weaker (Fig. 3A, B) as shown by Tomás et al. (2013) across a variety of genera. Contrarily, LMA was more strongly related to T leaf than to D leaf in Australian sclerophylls (Niinemets et al., 2009c). The median value of LMA for evergreen gymnosperms is higher (230 g m -2 ) than what was found here (214 g m -2 ), but the minimum and maximum values here fall into the gymnosperm range (Poorter et al., 2009) placing them at the low end of LES characterized by overall robust foliage and low photosynthetic capacity (Wright et al., 2004).
Prior this study the highest T cwm was recorded in ferns and allies: 0.81 µm (Tosens et al., 2016). The species here had even thicker mesophyll cell walls: up to 1.22 µm in P. nivalis. Both the average across the sample and across the gymnosperms (0.68 µm) were also very high compared with what has been previously found for non-stressed fully expanded leaves: from approximately 0.1-0.2 µm in annual to 0.4-0.5 µm in sclerophyllous species (Hanba et al., 2001;Terashima et al., 2006;Hassiotou et al., 2009;Tosens et al., 2012bTosens et al., , 2016Tomás et al., 2013). On the other hand, values of S c /S in our study were average in the majority of species as it varies from 1.6 to 32 m 2 m -2 worldwide (Terashima et al., 2006;Peguero-Pina et al., 2012;Tosens et al., 2016) and the lowest values have been recorded in ferns (in the majority of species S c /S<8) (Tosens et al., 2016), but across evergreen gymnosperms S c /S had only been measured in Abies alba and A. pinsapo (17 and 32 m 2 m -2 , respectively) (Peguero-Pina et al., 2012). However, CO 2 diffusion to Rubisco is greatly enhanced if the cell walls are fully lined with chloroplasts, with S c /S m =1, a condition that in practice is rare (Terashima et al., 2011). The ratio was remarkably high in the gymnosperms studied here, being over 0.96 for five species (see Supplementary Fig. S1A).
The LMA components thickness and density alter leaf photosynthetic capacity in opposite directions (Niinemets, 1999). Leaf thickness often reflects a greater number of chloroplasts, higher S c /S and S m /S for CO 2 diffusion, and greater concentration of photosynthetic machinery, while greater density is associated with lower net assimilation due to increased cell wall resistance (Hanba et al., 1999;Niinemets, 1999;Terashima et al., 2006). Contrarily, Slaton and Smith (2002) studied a wide range of genera, but did not find a correlation between S m /S and T mes . Similarly, no relationship was found between S c /S or S m /S and T mes across Australian sclerophylls (Tosens et al., 2012b). However, across the evolutionarily old species studied in this work, T mes correlated positively with S m /S, but regardless of the positive trend was not correlated with S c /S ( Supplementary Fig. S1B, C). However, the species with high S m /S also had higher S c /S ( Supplementary  Fig. S1A). This somewhat counterintuitive outcome results from the variability of S c /S due to the spore-bearing plants in   6. (A) The influence of chloroplast surface area exposed to intercellular airspaces (S c /S) on mesophyll conductance (g m/area ). (B-D) The influence of cell wall thickness (T cwm ) on (B) mesophyll conductance per mass (g m/mass ), (C) mesophyll conductance per area (g m/area ), and (D) the draw-down from intercellular airspaces to chloroplasts (C i -C c ). The data in (B) were fitted with a non-linear regression of the form y=-0.37ln(x)+0.09 and data in (C) of the form y=-26ln(x)+25. Data presentation and fitting are as in Fig. 3. the sample as their S c /S m ratio was lower than average. These relationships were not significant, when these two species were excluded from the analyses. S m /S and S c /S depended on mesophyll architecture rather than T mes or F ias due to variability in density, size and shape of mesophyll cells (e.g. lobes; Fig. 1C) at a given thickness as was previously suggested by Tosens et al. (2012a). Likewise, this resulted in LMA not correlating with S c /S (Fig. 3D).
Additionally, T cwm was not related to LMA (Fig. 3C) or density as is often assumed (e.g. Syvertsen et al., 1995) showing that T cwm may not vary with total leaf tissue density or other leaf tissue's cell wall thicknesses (Tosens et al., 2012b). Collectively, these results suggest that widely varying combinations of leaf anatomical traits occur at given values of LMA, and detailed anatomical studies are needed to relate g m to mesophyll structure.

Structural modifications of photosynthesis and mesophyll diffusion conductance
A area and g s were correlated and comparable to those found in spore-bearing plants and gymnosperms (James et al., 1994;Equiza et al., 2006;Franks, 2006;La Porta et al., 2006;Pardo et al., 2009;Zhang et al., 2015;Tosens et al., 2016). Flexas et al. (2012) set an average g m slightly above 100 mmol m -2 s -1 for conifers, and studies have shown a wide variation from as low as 18 up to 110 mmol m -2 s -1 in evergreen gymnosperms (Warren et al., 2003;De Lucia et al., 2003;Manter and Kerrigan, 2004;Black et al., 2005;Mullin et al., 2009;Peguero-Pina et al., 2012, 2016. Accordingly, g m varied from 10 to 124 mmol m -2 s -1 in this study. The values of C i -C c here fell in the higher part of the range found for angiosperms and were higher than previously found for conifers (Niinemets et al., 2005(Niinemets et al., , 2009bWarren, 2008;Hassiotou et al., 2009;Tosens et al., 2012b;Peguero-Pina et al., 2012). Importantly, g m influenced both C i -C c and photosynthesis significantly with stronger mass-based relations (Fig. 4).
Nitrogen content was similar to that previously found in conifers and positively related to A mass , as was also shown by Reich et al. (1995). M. glyptostroboides diverged in the analysis, which can be explained by its deciduous nature (Wright et al., 2002). A negative relationship between LMA and A mass has been shown previously for evergreens while area-based photosynthesis seems not to correlate with LMA, which is contrary to deciduous species, where higher LMA has been shown to correlate with higher A area (Reich et al., 1995;Wright et al., 2004). Accordingly, the same was found here (Fig. 5). Altogether, this illustrates strong investment in supportive structure (Niinemets et al., 2009c;Hassiotou et al., 2009).
Similarly to recent studies, the most important anatomical correlations with g m were S c /S and T cwm (Terashima et al., 2011;Tosens et al., 2012bTosens et al., , 2016Peguero-Pina et al., 2012;Tomás et al., 2013) (Fig. 6). The significance of each depends on foliage structure, e.g. in sclerophylls higher T cwm overshadows the influence of S c /S on g m (Tomás et al., 2013). The g m/mass -T cwm correlation was stronger than the area-based one (Fig. 6B, C) as g m/mass indicates the investment in cell walls and g m/area in the mesophyll area (Niinemets et al., 2009b).

The importance of structure as a limiting factor of photosynthesis
In order to overcome the various uncertainties related to g m estimations, it was additionally modelled from anatomy and estimated from A-C i curves (Ethier and Livingston, 2004). All three methods show similarly wide variation in g m across evolutionarily old species (Fig. 7). However, g m from anatomy tended to underestimate, whereas the A-C i curve method overestimated, g m for some species. When both methods were considered, g m from the variable J method was more strongly correlated with g m obtained from anatomy than from A-C i curves. Likewise, Tomás et al. (2013) found this across species exhibiting a wide variation of leaf morphologies. Yet, the model tends to over-or underestimate at the lower or higher end of g m (Tosens et al., 2012a,b;Peguero-Pina et al., 2012;Tomás et al., 2013;Fini et al., 2016). In this sample, g m from anatomy mostly underestimated for C. revoluta, C. sempervirens and P. abies (Fig. 7), while the model estimations where in strong agreement with g m from gas exchange for species with thick mesophyll cell walls and low g m . This is different from Tosens et al. (2016), where the model overestimated Fig. 7. Comparison of g m/area calculated based on Harley et al. (1992) with (A) g m/area calculated from anatomical measurements according to Niinemets and Reichstein (2003) and (B) g m/area calculated based on Ethier and Livingston (2004). Dashed lines represent 1:1 correlation. Data presentation and fitting are as in Fig. 3, but all species are included in the analyses. C. revoluta, C. sempervirens and P. abies were underestimated by the anatomical model.
for evolutionarily old Ophioglossum and Lycopodium species exhibiting very high T cwm (0.7-0.81 µm). The discrepancy between measured and modelled g m may arise from various uncertainties associated with both estimates. As for combined gas exchange-fluorescence estimation, g m is not the only component to affect C c : the amount of respiratory and   (Tholen et al., 2012). Concerning the predictive power of the anatomical model for species exhibiting a wide variation in T cwm , heterogeneous chloroplast thickness and morphology (Figs 1 and 2), the strongest uncertainties are related to the unknown variation of mesophyll cell wall porosity and the actual determinants of stroma resistance including the role of carbonic anhydrase (CA). Nobel (1991) proposed an upper value of cell wall porosity of 0.3 and that it varies linearly with T cwm . Indeed, 0.3 was used to calculate g m for Populus tremula with low T cwm and it gave good results (Tosens et al., 2012a). In a multispecies context the relationship between g m from anatomy and gas exchange was improved when porosity was modelled varying between 0.04 and 0.095 and being negatively correlated with T cwm (Tosens et al., 2012b;Peguero-Pina et al., 2012;Tomás et al., 2013). Some species here and in Tosens et al. (2016) exhibited very high T cwm compared with prior research, and therefore lower porosity values were used for species where T cwm >0.5 μm. However, as was emphasized by Tosens et al. (2012b), the linear decline of porosity with T cwm is hypothetical and the discrepancy between modelled and measured g m could be resolved by some other trait. Quantitative limitation analysis further confirmed that the mesophyll gas phase limitation to CO 2 diffusion is negligible as it accounted for less than 1% of the total limitation. This is even less than previously observed in a wide range of sclerophyllous angiosperm genera or in spore-bearing plants (Tosens et al., 2012a,b;Tomás et al., 2013). Using limitation analysis to further separate the contribution of the components of the liquid phase revealed that the limitation of T cwm and T chl had the greatest span, ranging from 6 to 72% and from 22 to 84%, respectively (Fig. 8B). This is in agreement with the negative correlation between g m/mass and T cwm (Fig. 6B). CO 2 faces the longest diffusion distance in the liquid phase after entering the chloroplasts . In this set of species, T chl varied between 1.97 and 5.36 μm (Fig. 2), which is thicker than previously reported explaining the high limitation by T chl . It is debated whether the CO 2 diffusion efficiency in the chloroplast stroma is mainly controlled by CA activity in the chloroplast  and if it varies across species (at normal cytosolic pH), being higher in evergreen sclerophylls in order to compensate for the low CO 2 diffusion efficiency through thick mesophyll cell walls (Gillon and Yakir, 2000), or whether CA is always sufficient and CO 2 interconversion is not a limiting factor of g m . Equally, this analysis adds to the evidence that CO 2 diffusion is importantly controlled by the physical diffusion distance within the chloroplasts (Tosens et al., 2012a,b;Peguero-Pina et al., 2012;Tomás et al., 2013).
Overall, the good fit between g m calculated from anatomy and gas exchange confirms that g m is explained to a large extent by inherent variations in mesophyll anatomy, while mesophyll cell walls and chloroplasts are the most important structural determinants of liquid phase conductance in species exhibiting high T cwm , T chl and heterogeneous chloroplast morphology. However, the limitation amplitude was large implicating versatility in its importance among different species (Fig. 8A). The smallest limitation on A area on average was biochemistry, which is similar to previous interspecies studies by Tomás et al. (2013) and Tosens et al. (2016). As the mesophyll and stomatal conductances were low, especially compared with maximal values measured so far (Wright et al., 2004;Flexas et al., 2012), it is expected that biochemistry plays a smaller role in limitations, because the photosynthetic enzymes are not saturated with CO 2 due to constraints on the previous parts of the pathway (Terashima et al., 2011;Niinemets et al., 2011). The limitations by stomata and mesophyll were similarly significant in the control of WUE i showing again that leaf mesophyll structure plays an important role in realized WUE i (Fig. 9).

The effect of divergence time on structure and physiology
The rate of photosynthesis depended on the age of the genera. This, however, was mostly due to the inclusion of spore-bearing plants reflecting the evolutionary increase in photosynthetic capacity from early plants to angiosperms (Brodribb et al., 2009;Carriquí et al., 2015;McAdam and Brodribb, 2015). However, when broad taxonomic groups are considered, species' evolutionary adaptation to light and water conditions can actually drive photosynthetic capacity more strongly than their evolutionary age (Tosens et al., 2016). In Fig. 9. Correlations of intrinsic water use efficiency (WUE i ) with (A) the share of photosynthesis limited by mesophyll conductance (l m ) and (B) share of photosynthesis limited by stomatal conductance (l s ). Data presentation and fitting are as in Fig. 3. this regard, it is important that A area depended on g m/area positively regardless of evolutionary age, showing that this correlation remained significant including the effect of divergence time in the model (see Supplementary Table S3a, b).
Analogously to A area , it has been suggested that there is an increase in g m/area through the evolution and diversification of embryophytes (Carriquí et al., 2015). While indeed gymnosperms studied here had a very low average g m/area , much lower than angiosperms with comparably tough foliage structure (e.g. Flexas et al., 2012), g m/area did not significantly depend on evolutionary age. As with within-group variability in ferns (Tosens et al., 2016) and angiosperms (Flexas et al., 2007b), this likely reflects individual species adaptation to specific habitat conditions.
Notwithstanding the lack of within-group evolutionary signal, low average g m/area and A area and the corresponding thick mesophyll cell walls for species with widely contrasting ecological strategies do support the preservation of old traits through evolution, suggesting apparent constraints on evolution. Different CO 2 /O 2 selection pressures at the time of divergence for ferns and angiosperms have been hypothesized to be the reason for lower average A area values for the former (Carriquí et al., 2015) as ferns evolved in several-fold higher CO 2 and slightly lower O 2 concentrations than angiosperms (Brodribb et al., 2009;Carriquí et al., 2015). In fact, in the high-CO 2 atmosphere where several of the thick-cell-walled species evolved (about 65-200 My ago), diffusional limitations exercised a lower control on the rate of photosynthesis. However, while this ancient trait has been preserved through evolution, g m has started to exercise an increasingly stronger control on foliage assimilation rates. Yet, this might be inevitable given the inefficient dynamic stomatal control under water-limited conditions potentially leading to excess plant water loss until hydraulic constraints force stomatal closure (Brodribb et al., 2005;Carnicer et al., 2013).

Conclusions
This study demonstrates important relationships between mesophyll structure and physiology in evolutionarily old plants. A large variability was found in ultrastructural traits determining LES relationships. Although LMA depended on leaf density and thickness, they were not correlated with T cwm or S c /S illustrating that LMA cannot be used as a universal explanation for photosynthetic structural constraints.
Mesophyll anatomy exerted major control on net assimilation rates. The three methods used in previous studies for interspecies comparisons were used to confirm g m values from gas exchange-fluorescence measurements and their anatomical nature. Although high variation was uncovered between species in g m , it was on average very low. This resulted mainly from the exceptionally thick cell walls, large chloroplasts, and low S c /S. The role of stroma deserves further attention.
These characteristics in evolutionarily older taxa support the preservation of ancient traits in evolution, although the role of functional adaptation and evolutionary constraint in leaf anatomy continues to be debated.

Supplementary data
Supplementary data are available at JXB online. Fig. S1. S c /S, S m /S and T mes correlations. Fig. S2. LMA and evolutionary age correlation. Table S1. Species' habitat and origin. Table S2. Dependence of WUE i on limitations. Table S3. A area and evolutionary age correlation.