Searching for a nearest living equivalent for Bennettitales: a promising extinct plant group for stomatal proxy reconstructions of Mesozoic pCO2

ABSTRACT To understand Earth´s climate variability and improve predictions of future climate change, studying past climates is an important avenue to explore. A previously published record of pCO2, across the Triassic–Jurassic boundary (TJB, ~201 Ma) of East Greenland, showed that Bennettitales (Anamozamites and Pterophyllum) responded in parallel to the empirically proven pCO2-responders Ginkgoales, reducing their stomatal densities by half across the TJB, indicating a transient doubling of pCO2. The abundance of fossil Bennettitales in Mesozoic strata and natural history museum collections worldwide offers enormous potential for further stomatal proxy pCO2 reconstructions, but a suitable nearest living equivalent (NLE) should ideally first be identified for this extinct plant group. Using specimens from herbarium collections, three species of cycads, historically considered the best NLE, were tested for pCO2 response, as well as two species of tree ferns, grown in experimental growth chambers. None responded to changes in pCO2, and were consequently rejected as NLEs. Finally, two species of ferns were selected from the literature, and produced very similar pCO2 compared to Ginkgoales. However, these understory ferns are not appropriate NLEs for Bennettitales due to differences in habitat and a distant evolutionary relationship. Future work should test additional plant groups, in particular seed plants such as basal angiosperms and Gnetales, for suitability as NLE for Bennettitales in pCO2 reconstructions, for example through biogeochemical fingerprinting using infrared microspectroscopy. Until an appropriate NLE is identified, Bennettitales pCO2 can be reconstructed based on cross-calibration of stomatal densities with those of co-occurring pCO2 responders, such as Ginkgoales.


Introduction
Carbon dioxide (CO 2 ) is a principal atmospheric greenhouse gas that has increased almost 50% in concentration since preindustrial times, and at an unprecedented rate during the last hundred years, due mainly to human activity. It is predicted to increase further and result in temperature increases of up to 4°C by the year 2100 (IPCC report, 2015), with severe socioeconomic and ecosystems impacts predicted. Reducing uncertainty in predictions of future climate change is one of today's greatest scientific challenges, with many important questions left unanswered. For example, though it is well established that CO 2 exerts radiative forcing that warms the planet, the exact relationship between atmospheric CO 2 concentration (pCO 2 ) and temperature, the so-called climate sensitivity, remains poorly constrained (Huber et al. 2014). Our principal tools for understanding the climate system, long-term pCO 2 records from geological archives and climate models, still lack the resolution to precisely constrain the role of pCO 2 in climate change, thus severely impeding our ability to project the consequences of future pCO 2 rise. Although marine proxy records of palaeo-temperature are considered highly reliable (e.g. Zachos et al. 2001;Pagani et al. 2011;Foster et al. 2012), marine pCO 2 records have recently come under some scrutiny for being strikingly variable and in some cases contradictory, or showing unexpected relationships with palaeo-temperature records, and this has ignited heightened interest in the terrestrial stomatal proxy for palaeo-pCO 2 reconstructions Beerling & Royer 2011;Barclay & Wing 2016;McElwain & Steinthorsdottir 2017).
Stomata are pores on plant leaf surfaces, through which carbon is acquired for photosynthesis from CO 2 and at the same time water vapour and oxygen are lost by diffusion. An inverse relationship exists between the density of stomata (SD = N stomata/mm 2 ), as well as the now more widely used stomatal index (SI = % stomata/all epidermal cells), and pCO 2 , caused by the physiological response of plants aiming to reduce the density of stomata when possible (when pCO 2 increases) in order to conserve water. This relationship was originally established based on herbarium material, showing that modern tree species have responded to the anthropogenic rise in pCO 2 by reducing SD significantly (Woodward 1987). The ratio between SD/SI in the past and at present is thus representative of the ratio between pCO 2 in the past and the known modern pCO 2 . Palaeo-pCO 2 is reconstructed from these interdependent ratios, either by direct calibration (the stomatal ratio method: McElwain & Chaloner 1995), or by constructing SD/SI response curves based on nearest living relative (NLR) leaf specimens collected from living plants or herbarium sheets, and/or from experiments, in all cases grown at known pCO 2 (transfer functions, now usually referred to as the stomatal index method, e.g. Barclay & Wing 2016). A more recent third methodology, gas-exchange modelling (Franks et al. 2014;Konrad et al. 2017), is sometimes characterised as more taxon-independent, but still must rely on measurements from modern plants that are closely related to the fossil plants in order to reconstruct pCO 2 . When no obvious NLR can be identified, such as for some extinct plant groups, palaeo-pCO 2 is usually reconstructed using the ratio between the SI of the fossil plants and the SI of a Nearest Living Equivalent (NLE), defined as a plant that is a morphological and/or ecological equivalent of the fossil plant, ideally a close or distant relative (McElwain & Chaloner 1995). Palaeo-pCO 2 is then inferred based on the known pCO 2 -SI relationship, using the SI ratio. The inverse relationship between SD/SI and pCO 2 has been repeatedly demonstrated for the majority of woody plant taxa, from disparate geological and ecological settings, from the Paleozoic (late Carboniferous) until today (see e.g. McElwain & Chaloner 1996;McElwain 1998;Barclay et al. 2010;Steinthorsdottir et al. 2011bSteinthorsdottir et al. , 2013Bai et al., 2015;Mays et al. 2015;Steinthorsdottir & Vajda 2015;Montañez et al. 2016;Steinthorsdottir et al. 2016aSteinthorsdottir et al. , 2016bMcElwain & Steinthorsdottir 2017;Wolfe et al. 2017;Slodownik et al. 2021;Steinthorsdottir et al. 2019aSteinthorsdottir et al. , 2019bSteinthorsdottir et al. 2020). The stomatal proxy is well supported by independent data and considered reliable in palaeo-climate reconstruction (see e.g. Beerling & Royer 2011;McElwain & Steinthorsdottir 2017).
The Bennettitales were a major group of extinct gymnosperms, which were widespread globally from the midlate Permian until their near-extinction at the end of the Cretaceous (with some perhaps surviving in refugia until the Oligocene, McLoughlin et al. 2010). Bennettitalean fossil leaves are very abundant worldwide in sedimentary deposits of Mesozoic age, frequently constituting a large proportion of the fossil flora, and are consequently well presented in many museum collections (see e.g. Harris 1932Harris , 1969Vakhrameev 1991;Watson & Sincock 1992;Knobloch & Kvacek 1997;Van Konijnenburg-van Cittert & Morgans 1999;Boyd 2000;Pott et al. 2007;Crane & Herendeen 2009;Pott & McLoughlin 2009;McLoughlin et al. 2010;Pott & McLoughlin 2014). This important group thus represents some of the most common Mesozoic fossil plants, and >900 Bennettitales macro leaf specimens are hosted in the collections at the Swedish Museum of Natural History, principally from Skåne (Scania) in southern Sweden as well as from Greenland and Australia (Fig. 1A-F).
Bennettitalean leaf cuticles have previously been recognized as reliable proxies for palaeo-pCO 2 concentrations, based on their parallel SI response to those of fossil Ginkgoales across the Triassic-Jurassic boundary (TJB) of Sweden, Greenland and the UK (McElwain et al. 1999;Steinthorsdottir et al. 2011aSteinthorsdottir et al. , 2011b. The Ginkgoales are a group of empirically well-established pCO 2 -responders, including the single extant species, the "living fossil" Ginkgo biloba (e.g. Beerling et al. 1998;Chen et al. 2001;Haworth et al. 2015;Barclay & Wing 2016;Retallack & Conde 2020), and Ginkgoales fossil leaves have been used extensively for Mesozoic and Cenozoic pCO 2 reconstructions (Chen et al. 2001;Royer 2003;Sun et al. 2007;Quan et al. 2009;Smith et al. 2010;Barbacka 2011;Steinthorsdottir et al. 2011b;Mays et al. 2015;Sun et al. 2016;Wu et al. 2016;Sun et al. 2018;Zhou et al. 2020). Although fossil Ginkgoales, having a nearest living relative, may be considered superior in palaeo-pCO 2 reconstructions, they are not always abundantly present in museum collections or from geological strata of interest, and the potential to use alternative fossil plants, such as the Bennettitales, may enable more numerous and/or higher-resolution pCO 2 reconstructions across important climate change transitions in the past. In addition, Bennettitales (as well as Ginkgoales) cuticle morphology may indicate volcanic pollution (primarily by SO 2 ) in the past, thus recording important environmental signals linked to climate and environmental change as well as mass extinctions (Haworth et al. 2012;Steinthorsdottir et al. 2018).
In order to fully harness the potential of Bennettitales as pCO 2 recorders for the Mesozoic using the stomatal proxy method, it would be useful to identify an appropriate NLE. In this context, a reliable NLE would show a significant response to changes in pCO 2 , in addition to being morphologically and ecologically comparable to the extinct proxy group. Further, to mask the effects of species-specific stomatal responses, including several suitable NLE candidates and using the combined stomatal index may be considered the superior approach (e.g. Steinthorsdottir et al. 2019aSteinthorsdottir et al. , 2019bSteinthorsdottir et al. 2020). Historically Bennettitales, also referred to as the Cycadeoideales in some classifications, have been compared to Cycadophyta (cycads), based on similar macro-leaf morphology and co-occurrence in the fossil record. The two have even been referred to collectively as the Cycadophytes (Norstog & Nicholls 1997). All proxy methods that use morphological traits when direct measurements are not available assume strong physiological conservatism between extant and fossil groups (because they look the same, they work in the same way). However, this assumption is usually not tested, which may weaken the proxy approach. Here, the potential as NLEs for Bennettitales of three cycad species (Lepidozamia peroffskyana, L. hopei and Zamia furfuracaea), and two tree fern species (Dicksonia antarctica and Cyathea cooperi), all selected based on their ecological and/or morphological similarity to the East Greenland Bennettitales used in the published pCO 2 reconstruction of Steinthorsdottir et al. (2011aSteinthorsdottir et al. ( , 2011b, was tested.

Selecting potential NLEs to Bennettitales
As outlined above, it was established as a criterion that any NLE to Bennettitales should be a pCO 2 -responder, reducing stomatal frequency with rising pCO 2 . Herbarium material from three species of cycads was recorded in the first trial, and living material from two species of tree ferns, grown in experimental chambers under ambient and elevated pCO 2 treatments, in the second. All microscopy slides and experimental leaf material reported in this study are hosted at the Plant Palaeoecology and Palaeobiology Research group, Department of Botany, Trinity College, Dublin, and are available for further study on request. The herbaria collections specimen numbers, when assigned, are listed in the supplementary information (SI Table 1).

Cycads
Cycads have until recently been considered the most appropriate NLEs of Bennettitales due to their macromorphological similarity and were, therefore, chosen for testing in the first instance. In particular, Lepidozamia peroffskyana Regel and L. hopei Regel ( Fig. 2A-B) were considered good ecological equivalents, due to their habitat preferences similar to that of the Bennettitales from the Triassic-Jurassic of East Greenland (understory to mid-canopy in wet environments, often river flood plains) and morphological similarities (leaves divided into leaflets without midrib) (Jones 2002;Falcon-Lang & Cantrill 2002;McElwain et al. 2007). Both species of Lepidozamia are Southern Hemisphere cycads: this is also in accordance with the East Greenland palaeo-plant community, which has been inferred to most closely resemble a Southern Hemisphere broad-leaf conifer forest (McElwain et al. 2007). Zamia furfuracea L.f. (a Mexican cycad adapted to drought and thus not accurately reflecting the palaeoecological conditions in which the Bennettitales grew) was also included in the group of NLEs as a good morphological NLE. This is due to the leaf morphology, which most closely reflects bennettitalean foliage architecture (here Anamozomites and Pterophyllum,   1A-B, F), with its shorter and wider leaflets than most other cycads. Together, the three cycad species were thus hypothesised to be the best possible candidates as NLEs of Bennettitales across the TJB in East Greenland.

Tree ferns
Two species of tree ferns, Dicksonia antarctica and Cyathea cooperi (Fig. 2C-D), were additionally tested as suitable NLEs of Bennettitales. Although not as readily associated with Bennettitales as the cycads, these tree ferns share many traits with whole-plant reconstructions of some Bennettitales, such as a short trunk, pinnate leaves and an ecological position as mid-canopy to understory plants (Falcon-Lang & Cantrill 2002;Hunt et al. 2002;McElwain et al. 2007). Ferns are also ancient plants ("living fossils"), with a long fossil record (Schneider et al. 2004). In addition, although cycads may be more similar to Bennettitales macromorphologically, ferns and Bennettitales are more similar micromorphologically, both having cuticles with crenulated epidermal cell walls and similar stomatal appearance (Figs. 1C, 2C-D). A leaf from each of three plants of D. antarctica and C. cooperi respectively was collected from plants grown in either ambient (380 ppm) or elevated (1500 ppm) pCO 2 as a pilot study to establish whether or not these species respond to changes in pCO 2 .

Leaf database
A database of 43 cycad leaflets, comprising 24 leaflets from Lepidozamia peroffskyana herbarium specimens, ten from L. hopei and nine from Zamia furfuracea (see SI Table 1: All data) were obtained by sub-sampling herbarium sheets at the Royal Botanic Gardens Kew Herbarium, London, UK (Kew) and requesting loans of additional sub-sampled specimens from the Australian National Herbarium, Canberra, Australia (CANB) and the Field Museum of Natural History, Chicago, USA (FMNH). One published Z. furfuracea SI value was added to the analysis dataset (McElwain et al. 1999). The specimens range from the 1780's to the late 1990's and thus span an anthropogenic rise in atmospheric pCO 2 of >30% or ~90 ppm (from ~280 ppm to ~370 ppm (Keeling et al. 2005; CO 2 data from noaa.org (1800-1850 and 1958-present) and data.giss. nasa.org ). To compare, we collected six leaves from each of the tree ferns (Dicksonia antarctica and Cyathea cooperi), which were grown in the UCD PÉAC facility's controlled atmosphere walk-in growth chambers (CONVIRON BDW40) in pCO 2 concentrations of 380 ppm and 1500 ppm, to assess their responsiveness to pCO 2 (SI Table 1).

Laboratory and statistical methods
For the cycad leaflet database, circular discs of approximately 1 cm in diameter were extracted from the mid lateral section of each herbarium specimen leaflet using a hole-puncher. The leaf discs were then macerated using 50/50 glacial acetic acid (CH 3 COOH) and hydrogen peroxide (H 2 O 2 ), following the method of . When all mesophyll was dissolved, the upper and lower cuticles were separated in a Petri dish with water using small brushes under the dissection microscope. The lower (abaxial -stomata bearing) cuticle discs were stained with a staining agent (Sudan IV) and mounted on permanent slides using Keiser's gel (glycerin gelatin). For the tree fern leaf database, the cuticular morphology was recorded using the dental impression/nail polish technique, which is extremely fast and safe. Roughly equal amounts of dental base and catalyst were first mixed together (both Coltène® "President" polyvinylsiloxane), this combination dried to a rubbery solid substance, replicating leaf surface morphology in fine detail. Before hardening, each leaf was pressed about halfway down into the mixture and left overnight, when the leaves were carefully removed. Next, 2-3 coats of lightly-tinted clear nail polish were applied onto the casts made by each leaf. In approximately 20-30 minutes, when the nail polish had dried but not completely hardened, the coat of nail polish was removed and placed cast-side up on a glass slide beneath a cover slip. Each macerated cuticle disc and nail polish slip was photographed 7-10 times at different sites at ×200 magnification using a Leica camera (Leica DFC300 FX), mounted on a light microscope (Leica DM 2500). Auto-Montage Pro (©Synoptics, version 5.03.0061), a software that stacks images of specimens possessing relief in the z plane, was used to create a single flat image when necessary. The images were engraved with 0.09 mm 2 grids using the software AcQuis (©Synoptics, version 4.01.10), stomata and epidermal cells were then counted within each grid using the software ImageJ (1.39 u, NIH, USA), following the method of Poole and Kürschner (1999), and the cumulative mean SI for each counted specimen was recorded. SI data was entered into Excel spreadsheets and the standard error calculated (see SI Table 1). For the cycad analyses, we employed linear regression, with the SI as dependent variable and the ppm of CO 2 in the year of collection as explanatory variable. Three independent models were fitted for the three species. For the tree fern experiment, we analyzed the data using Analysis of Variance. A model with SI as dependent variable and species, pCO 2 treatment (elevated or ambient), and the interaction between species and pCO 2 treatment was used. In both analyses, assumptions of the models (normality of the residuals, homoscedasticity, absence of influential points) were checked by eye by plotting theoretical quantiles vs the standardized residuals (normality of the residuals), fitted values vs the square root of the standardized residuals (heteroscedasticity), and leverage vs the standardized residuals (influential points). All statistical analyses were conducted in R (R Core Team 2019).

Cycads: herbarium study
The tested cycad species did not respond by decreasing their SI at increased pCO 2 . Lepidozamia peroffskyana showed no change in their average SI of 3.8% across the ~90 ppm pCO 2 increase (p-value 0.7495, R 2 0.00473) (Fig. 3). Lepidozamia hopei also showed no statistically significant change in their average stomatal index of 3.6% (p-value 0.885, R 2 0.0028). Zamia furfuracea showed a statistically significant positive correlation to rising pCO 2 , from an SI of approximately 7.5% at preindustrial atmospheric pCO 2 to approximately 10% in 1998 at ~367 ppm pCO 2 (Fig. 3 p-value 0.0002, R 2 0.8354), contrary to the response assumptions of the stomatal proxy theory.

Tree ferns: growth chamber study
The two tree ferns species, Dicksonia antarctica and Cyathea cooperi ( Fig. 2E-F), were grown in the PÉAC controlled atmosphere chambers at 380 ppm (ambient pCO 2 treatment) and at 1500 ppm (elevated pCO 2 treatment). Three plants per species per treatment were grown from seedlings, and the leaves collected were fully developed in their respective atmospheric concentration (as per Lake et al. 2001). Our ANOVA showed no influence of treatment on SI (F-value 2.171, p-value 0.17882) or the interaction between species and treatment (F-value 0.845, p-value 0.3848), while it found significant difference between the SI of the two species (F-value 19.246, p-value 0.0023) (Fig. 4). The two tree fern species were consequently found not to be suitable NLEs to Bennettitales, since they do not reduce their SI with increasing pCO 2 .

Ferns: literature study
In the search for applicable NLEs for Bennettitales for the original pCO 2 reconstruction (Steinthorsdottir et al. 2011b), we looked for additional potential NLEs outside of the cycad and tree fern studies conducted. A literature search unveiled a promisingly appropriate SI of ~24% (reflecting the higher SI values of Bennettitales versus Ginkgoales in TJB material) by combining two species of ferns that had been suggested as potential pCO 2 responder NLEs; the proven pCO 2 -responder Stenochlaena palustris (Beerling et al. 2002) and the pseudo tree fern Todea barbara. The combined SI of the two NLEs selected here produces TJB pCO 2 values remarkably similar to those calibrated using Ginkgo biloba as NLE for Ginkgoales in the same section. Stenochlaena palustris recorded SI of 31.9% at 360 ppm (Beerling et al. 2002), whereas T. barbara recorded SI of 16.2% at modern CO 2 (here ~390 ppm), resulting in a combined NLE SI of 24% at ~375 ppm. Stenochlaena palustris (Blechnaceae) is an understory fern with leaves that are only once-dissected, displaying a leaf morphology somewhat similar to Bennettitales. It is a primary colonizer in warm and humid Indo-Malaysian and African vegetation, sometimes as a climbing fern (Wolfe & Upchurch 1987). Stenochlaena palustris is a proven pCO 2 -responder (Beerling et al. 2002) and in combination with ecological and morphological considerations partially satisfies the criteria for a Bennettitales NLE. Todea barbara (Osmundaceae), although not a true tree fern, has a short trunk and thus occupies a low mid-canopy position. It is a Southern Hemisphere fern, often found in moist, riverside habitats (Page 1979). The arborescent habit of these ferns may make them potential NLE for Bennettitales, since many bennettitalean plants, e.g. Wielandiella, Williamsoniella, and Becklesia were arborescent (Watson & Sincock 1992;Pott & McLoughlin 2014). In addition, ancestral Osmundaceous ferns, such as Cladophlebis and Todites, have been found in East Greenland (and other equivalent geological sections), in comparable vegetational communities to the Bennettitales studied here (M. Popa, pers. comm.), indicating similar ecological preferences. Although not tested, T. barbara is tentatively hypothesized to be a pCO 2 responder, since Osmundaceae includes the empirically proven pCO 2 -responder Osmunda regalis (Wagner et al. 2007;Haworth et al. 2013). Todea barbara is morphologically (whole-plant) and ecologically a reasonably good analogue to Bennettitales, and the presence of the phylogenetically related Osmundaceae foliage in association with Bennettitales is an additional argument for its, at least partial, suitability as NLE. However, despite some similarities as outlined above and good SI fit for the TJB pCO 2 reconstruction of Steinthorsdottir et al. (2011b), ferns are only distantly related to Bennettitales and more appropriate NLEs remain to be identified.

Discussion
The first target NLEs for Bennettitales, the three species of cycads from herbarium collections, did not show the expected response of decreasing their SI with increasing pCO 2 . These results are independently supported by a separate study on living L. peroffskyana grown in ambient (380 ppm) and elevated (1500 ppm) pCO 2 , recording no SI response for this species (Haworth et al. 2013). In that study, SI for L. peroffskyana grown for eighteen months at 380 ppm pCO 2 was 4.6%, compared to 4.8% for cycads grown at 1500 ppm pCO 2 (Haworth et al. 2013). A significant increase in L. peroffskyana stomatal size (pore length) of 20.4% was recorded, however, under elevated pCO 2 relative to ambient (Haworth et al. 2013), and this may indicate an alternative response of cycads to pCO 2 . This, again, does not resemble the response of the TJB Bennettitales -their stomatal pore length instead decreased significantly when reconstructed pCO 2 was most elevated, to ~2000 ppm (Steinthorsdottir et al. 2012). The findings indicate that cycads are non-responders in terms of altering stomatal numbers to increases in atmospheric pCO 2 in both natural and experimental conditions, and should not be used for Bennettitales-based palaeo-pCO 2 reconstructions. The low SI-numbers further suggest fundamental differences between the pCO 2 responses of Bennettitales and cycads, since fossil plants have repeatedly been shown to have lower SI in the past (pre-Paleogene) than their relatives do now (e.g. McElwain 1998;Retallack 2001), whereas Bennettitales have significantly higher pre-Paleogene SI than cycads do. The cycads tested here are, therefore, not suitable as NLEs for Bennettitales.
Cycads were long thought to be good or relatively good NLE for Bennettitales, due principally to their co-occurrence in Mesozoic assemblages, the assumed long-lineages of both groups, as well as the morphological similarities between their leaf fossils, making them the original primary target when searching for NLEs for Bennettitales pCO 2 reconstruction using the stomatal ratio method. Cycadales was until recently considered an ancient plant group (Jones 2002), with preserved biological and ecological function (Haworth et al. 2011). Modern cycads have historically been referred to as "living fossils" due to their similarity to their fossil ancestors, and morphological continuity has typically been associated with physiological conservatism in palaeobotanical research. In the case of cycads, however, a series of more recent studies indicate that extant cycads are in fact considerably more derived than previously assumed, and may therefore not be reliable representatives of Mesozoic Cycadales, much less Bennettitales. Molecular dating analyses found that most of the modern species diversity of the Cycadales originated in the late Miocene (Nagalingum et al. 2011). Many fossil taxa that were originally identified as close relatives of extant genera have been identified as potential members of entirely extinct lineages, with little shared evolutionary history with the modern species (Erdei et al. 2019;Barone Lumaga et al. 2015). The few fossil cycads for which a relationship between fossil and extant taxa has been validated by phylogenetic analyses (Coiro & Pott 2017) show no evidence of a drastic response of SD or stomatal size to pCO 2 (Hill et al. 2019). The different response of cycads and Bennettitales to pCO 2 is thus not surprising in the light of morphological phylogenetic analyses showing that the two groups are only distantly related (Crane 1985;Doyle & Donoghue 1986;Doyle, 2008;Hilton & Bateman, 2006;Coiro et al. 2018). These results have been also validated by Fouriertransform infrared spectroscopy (FTIR) of leaf cuticles (D'Angelo & Zodrow 2015;D'Angelo et al. 2010;Vajda et al. 2017). Additionally, a recent study by Jardine et al. (2019) examined for the first time the impact of pCO 2 on Ginkgo cuticle chemistry, an avenue that should be further explored with Bennettitales cuticles in the future. Morphological phylogenetic analyses favour a close relationship between Bennettitales and angiosperms with or without Gnetales (Friis et al. 2007;Doyle, 2008;Coiro et al. 2018 Bennettitales stomata responded in parallel with those of Ginkgoales across the TJB by reducing their SI by halfa response that necessitates a comparison to an NLE to semiquantify pCO 2 . The Triassic-Jurassic Ginkgoales had an obvious NLE in the form of the only surviving species in the clade, Ginkgo biloba, another "living fossil" and a proven pCO 2 responder. Here we tested cycads and tree ferns as potential NLEs for Bennettititales, but none of the tested species responded to elevated pCO 2 by reducing their SI. The underlying principle of the stomatal proxy method is that plants seek to improve their water use efficiency (WUE), by minimizing water vapor-loss through reduction of stomatal densities, in times of abundantly available CO 2 in the atmosphere (high pCO 2 ), a response that a large majority of woody plants employ (Woodward 1987;Royer 2001;McElwain & Chaloner 1995;McElwain & Steinthorsdottir 2017).
The fact that the three tested species of cycads and two of tree ferns did not respond in the predicted manner to changes in pCO 2 indicates that they either employ a different mechanism to the majority of plants (including the extinct Bennettitales) to secure maximum WUE, for example by altering stomatal conductance, and/or that they have not evolved to maximize WUE, due to water being a non-limited resource. For example, the cycad Dioon sonorense shows high levels of WUE in the field, thus suggesting the presence of alternative mechanisms in this group (Álvarez-Yépiz et al. 2017). Furthermore, a study on the response of L. peroffskyana to increasing pCO 2 revealed that it does indeed respond but not by altering stomatal numbers; instead, it reduces stomatal conductance (G s ; mmol m −2 s −1 ) (Haworth et al. 2013). As pCO 2 was increased from 400 ppm to 2000 ppm, plants responded rapidly with pronounced reductions in G s (−58.5%), indicating that this species employs immediate and active stomatal control in response to environmental signals such as pCO 2 , in contrast to other species, including Bennettitales and Ginkgoales, that employ the slower method of amending stomatal densities. This rapid response was further confirmed in a later study on the rate of stomatal closure in response to darkness in an evolutionary range of plants from ferns to angiosperms (Elliott-Kingston et al. 2016). It has generally been assumed that small stomata respond faster to environmental signals than large stomata and the study was designed to test this assumption. Of the seven species tested, the species with the largest stomata, L. peroffskyana (35.6 µm), had the second fastest closing time in response to darkness (after Hordeum vulgare, a graminaceous angiosperm with dumbbell-shaped stomata), confirming that this species responds to environmental triggers through active physiological control of stomatal opening and closing, rather than through altering stomatal density (Elliott-Kingston et al. 2016). The discrepancy in physiological response of these extant plants and the extinct Bennettitales to rising pCO 2 indicates that neither the cycads L. peroffskyana, L. hopei and Z. furfuracea, nor the tree ferns D. antarctica and C. cooperi, should be employed as NLEs for Bennettitaleans when reconstructing palaeoatmospheric-pCO 2 . Furthermore, even though the fern (here S. palustris and T. barbara) SI data employed from the literature did result in comparable TJB Bennettitales-based pCO 2 to Ginkgoales-based pCO 2 , the distant evolutionary relationship of these seedless plants to Bennettitales and related groups suggests that, despite some macro-and micro-morphological similarities, these too should be discounted as appropriate NLEs.
Our results therefore indicate that other plant groups should be explored in future studies as better NLEs for the Bennettitales. The prospect of employing Gnetales as NLEs for the Bennettitales is intriguing. However, the morphology and ecology of the two genera of the Gnetales which share paracytic mesogenous stomata with the Bennettitales (Rudall & Bateman 2019; Rudall & Rice 2019) might be too divergent. Welwitschia, a plant growing in the Namib desert, has continuously growing linear leaves, while Gnetum has angiosperm-like leaves with multiple vein orders, and inhabits the understory of tropical forests. Early angiosperms, such as ANA-grade angiosperms, Chloranthaceae and mangoliids, might offer some insights. However, the diversity in epidermal anatomy and ecology in these groups would probably require a comparative phylogenetic framework to disentangle the different responses and find an appropriate NLE for Bennettitales.
The possibility also remains that even if the group most closely related to Bennettitales is firmly identified, it may not be an appropriate NLE for Bennettitales stomatal proxy pCO 2 reconstruction. The vast museum collections and abundance in geological strata of Mesozoic Bennettitales nonetheless provide an important opportunity to explore Earth system climate sensitivity by reconstructing palaeo-pCO 2 , for example as relative change recorded at high resolution across important climate change episodes. Another way forward is to compare the stomatal densities of Bennettitales to those of the proven pCO 2 -responders Ginkgoales when co-occurring in multiple coeval stratigraphical levels, and to cross-calibrate using the differences in SD/SI between the two taxa to reconstruct pCO 2 . This methodology has previously been successfully employed to reconstruct Late Triassic pCO 2 using fossil leaves of the extinct seed fern Lepidopteris ottonis from Germany (Bonis et al. 2011), as well as from a database of museum specimens from Scania, Sweden, hosted at the Swedish Museum of Natural History (Slodownik et al. 2021). The cross-calibration approach provides a viable alternative to NLE-based stomatal ratio palaeo-pCO 2 reconstruction, and should be explored in future work.

Conclusions
Reconstructing pCO 2 during past climate change episodes is an important tool to understand the workings of the Earth system and better predict the path and consequences of future anthropogenic climate change. The stomatal densities of Bennettitales responded in parallel to those of proven pCO 2responder Ginkgoales to a transient doubling of pCO 2 across the Triassic-Jurassic boundary of East Greenland. The potential of Bennettitales fossil leaves in palaeo-pCO 2 reconstruction is of considerable interest, given the abundance of these fossils in Mesozoic strata and in museum collections worldwide. To calibrate pCO 2 using stomatal densities in the stomatal proxy method, an NLE must ideally be selected for the fossil plants. An appropriate NLE should be a pCO 2 responder with comparable morphological and/or ecological characters to the fossil plants investigated. Here, three species of cycads, the group historically considered morphologically and ecologically closest to Bennettitales, were first tested as potential NLEs, using herbarium material spanning the recent anthropogenic rise in pCO 2 , but were found to be unresponsive or weakly positively correlated to pCO 2 , and were rejected. Two species of tree ferns were tested next, of ancient groups somewhat morphologically similar to Bennettitales, grown in experimental chambers under ambient and elevated pCO 2 , but were also unresponsive, and consequently rejected. Reconstructed pCO 2 using published SI data of two additional fern species as Bennettitales NLEs produced similar results to Ginkgoales pCO 2 in a previous study, but are not considered appropriate NLEs either, due to distant evolutionary relationship and ecology. Future work may explore the potential of Gnetales and/or basal angiosperms as Bennettitales NLEs. Bennettitales fossil leaves may still be utilized for high-resolution stomatal proxy pCO 2 reconstructions, by cross-calibration with other cooccuring, coeval pCO 2 responders with known NLEs.