Rapid response to anthropogenic climate change by Thuja occidentalis: implications for past climate reconstructions and future climate predictions

Carbon isotope values of leaves (δ13Cleaf) from meta-analyses and growth chamber studies of C3 plants have been used to propose generalized relationships between δ13Cleaf and climate variables such as mean annual precipitation (MAP), atmospheric concentration of carbon dioxide ([CO2]), and other climate variables. These generalized relationships are frequently applied to the fossil record to create paleoclimate reconstructions. Although plant evolution influences biochemistry and response to environmental stress, few studies have assessed species-specific carbon assimilation as it relates to climate outside of a laboratory. We measured δ13Cleaf values and C:N ratios of a wide-ranging evergreen conifer with a long fossil record, Thuja occidentalis (Cupressaceae) collected 1804–2017, in order to maximize potential paleo-applications of our focal species. This high-resolution record represents a natural experiment from pre-Industrial to Industrial times, which spans a range of geologically meaningful [CO2] and δ13Catm values. Δleaf values (carbon isotope discrimination between δ13Catm and δ13Cleaf) remain constant across climate conditions, indicating limited response to environmental stress. Only δ13Cleaf and δ13Catm values showed a strong relationship (linear), thus, δ13Cleaf is an excellent record of carbon isotopic changes in the atmosphere during Industrialization. In contrast with previous free-air concentration enrichment experiments, no relationship was found between C:N ratios and increasing [CO2]. Simultaneously static C:N ratios and Δleaf in light of increasing CO2 highlights plants’ inability to match rapid climate change with increased carbon assimilation as previously expected; Δleaf values are not reliable tools to reconstruct MAP and [CO2], and δ13Cleaf values only decrease with [CO2] in line with atmospheric carbon isotope changes.


INTRODUCTION
The concentration ([CO 2 ]) and isotopic value (d 13 C atm ) of atmospheric CO 2 are changing at a pace unprecedented in geologic time (Keeling et al., 2005;Zhang et al., 2013).
These changes have been accompanied by regional changes in mean annual temperature (MAT), mean annual precipitation (MAP), maximum summer temperature, and other climate variables (Yonetani & Gordon, 2001). The rapid decline in the carbon isotopic composition of CO 2 (d 13 C atm ) due to fossil fuel combustion, deforestation, and other human inputs, is known as the Suess Effect, and is a chemical representation of anthropogenic changes to the atmosphere-and more broadly, the environment. d 13 C atm values provide a useful way to see changes in CO 2 sources, sinks, and fluxes in the modern environment (Keeling, 1979;Boutton, 1991;Deines, 1992). It can also be applied to geologic problems (Schmitt et al., 2012) due to the naturally differing isotopic compositions of different CO 2 sources (e.g., methane, volcanism). d 13 C atm values are particularly useful because they are parameters in models that reconstruct past changes to atmospheric [CO 2 ] using paleosol carbonates (Cerling et al., 1991(Cerling et al., , 1992 or atmospheric [CO 2 ] using plant stomatal parameters (Franks et al., 2014). Direct measurements of d 13 C atm values only go back 50 years due to technological limitations, and longer-reaching ice core CO 2 bubbles (∼800,000 years) are poorly resolved for recent times and limited by the presence of glacial ice (Keeling & Whorf, 2004;Augustin et al., 2004;Barnola et al., 1987;Trudinger et al., 1999;Petit et al., 1999). The biosphere provides an excellent system that directly interacts with the atmosphere and fills the gap to provide high-resolution recent and long-term records, potentially extending into geologic time (Arens, Jahren & Amundson, 2000).
This direct interaction means that plants potentially provide a robust record of d 13 C atm values in their own leaf carbon isotope values (d 13 C leaf ) and fractionation values (D leaf , Eq. (1); Farquhar, Ehleringer & Hubick, 1989;Feng, 1999;Farquhar & Sharkey, 1982), which gives insight into changes in carbon assimilation over time. In Eq. (1), a represents the fractionation of d 13 C due to diffusion in air (4.4‰) and b represents the fractionation due to the carboxylation (instigated by the Rubisco enzyme, 27‰; Farquhar, Ehleringer & Hubick, 1989). These fractionation factors are compiled and multiplied by the ratio of C i (intercellular [CO 2 ]) to C a (atmospheric [CO 2 ]), a ratio that is often used to represent water use efficiency.
While a and b are thought to be constant, we know that d 13 C atm and C a are changing rapidly. This could result in a corresponding change in D leaf values as plants adapt to increased [CO 2 ] or subsequent regional climate changes, for example, systematic changes in local precipitation). Alternatively, D leaf values of leaves may stay constant but show marked changes in d 13 C leaf values corresponding to changes in d 13 C atm values if leaves are incorporating d 13 C atm into leaf tissues at a rate unaffected by other climate conditions.

CARBON ISOTOPES RELATED TO CLIMATE VARIABLES
Previous studies have related D leaf values to climate variables such as MAP, water availability and soil moisture (Diefendorf et al., 2010;Kohn, 2010;Wernerehl & Givnish, 2015;Mårtensson et al., 2017), MAT (Troughton & Card, 1975;O'Leary, 1993); latitude (Diefendorf et al., 2010;Kohn, 2010), [CO 2 ] (Schubert & Jahren, 2012, altitude (Korner, Farquhar & Wong, 1991;O'Leary, 1993), seasonality (Ehleringer, Phillips & Comstock, 1992), and d 13 C atm values during growth seasons (Peñuelas & Azcón-Bieto, 1992;Arens, Jahren & Amundson, 2000;Pedicino et al., 2002). The studies that incorporate potential influence from a wide range of climate variables have been conducted via meta-analysis with no normalized collection procedure or investigated species, or via growth chamber experiment conducted under idealized conditions. The few studies that have used naturally-obtained specimens (i.e., natural history collections such as herbaria) to look at isotope change over time and changing atmospheric conditions ([CO 2 ] and d 13 C atm values) have focused on localized regions with little range in climate (i.e., all dry, mid-to high-altitude, hot regions of eastern Arizona/western New Mexico, or the Mediterranean climate of Catalonia; Peñuelas & Azcón-Bieto, 1992;Pedicino et al., 2002). These collections-based experiments are limited in scope and while they provide information on specific ecosystems, do not address these biosphere-atmosphere interactions across climate regimes or on a regional and global scale. Very little is known about whether individual species respond to any, some, or all of these potential forcings across a range of climatic conditions.

POTENTIAL CLIMATE DRIVERS [CO 2 ] and elevation
We would expect higher [CO 2 ] to affect biochemical discrimination because of the known effects elevated [CO 2 ] has on stomata (size, density, and conductance; Woodward, 1987;Woodward & Bazzaz, 1988;Tognetti et al., 2000;Ainsworth & Rogers, 2007). In a meta-analysis of trees in European temperate and boreal forests, leaves responded to an increase in [CO 2 ] with a significant (21%) decrease in stomatal conductance (the rate of passage of atmospheric CO 2 into plant tissue; Medlyn et al., 2001). Increased [CO 2 ] also causes a decrease in stomatal density (the number of pores on a leaf surface) and stomatal index (the number of pores compared to the number of total epidermal cells); these indices have been liberally applied to the geologic record to reconstruct [CO 2 ] (Retallack, 2001;Beerling & Royer, 2002;Roth-Nebelsick et al., 2014). Because stomata directly control the flow of carbon dioxide into leaves and control carbon isotope fractionation by diffusion (Farquhar, Ehleringer & Hubick, 1989), changes in stomatal parameters could affect fractionation as well. Additional factors must be considered; the dependence of D leaf on [CO 2 ] may in part be due to isotopic discrimination associated with photorespiration (Schubert & Jahren, 2018). Indeed, previous growth chamber studies in prescribed CO 2 environments showed increased carbon isotope fractionation with increased [CO 2 ]; Schubert & Jahren (2012) found a strong hyperbolic correlation (r > 0.94) between [CO 2 ] and D leaf values in two species of herbaceous angiosperms. Based upon these growth chamber experiments (Schubert & Jahren, 2012, the relationship between [CO 2 ] and D leaf is expected to be most sensitive at geologically low [CO 2 ] (including pre-Industrial to present values) as it was in the levels present during plant growth in this study.
Elevation has been shown to factor into carbon isotope discrimination but is frequently not evaluated independently, due to its covariant relationship with climate variables such as temperature, vapor pressure, partial pressure of CO 2 (pCO 2 ), soil [CO 2 ], and soil texture (Diefendorf et al., 2010). A study looking at Salix herbacea leaves along an altitudinal gradient (2,000-2,800 m) in Austria showed a decrease in carbon isotope value with increased altitude, did not account for corresponding changes in other climate variables (Beerling, Mattey & Chaloner, 1993). Another study done in Utah and New Mexico using a number of desert and woodland species, including angiosperms and gymnosperms, found similar negative trends in d 13 C leaf with increased altitude without controlling for other climate variables (Van de Water, . In 2010, a meta-analysis assessing carbon isotope fractionation and discrimination values across a wide range of C 3 plants (Diefendorf et al., 2010) found that when combined with MAP, elevation explained 61% of variability in D leaf values. Based upon these studies, we included elevation as a potential variable in our study, however, we chose sample locations to minimize changes in elevation because it is difficult to completely separate this variable from regional variation in [CO 2 ] and d 13 C atm . δ 13 C atm With the current increase in [CO 2 ], we have observed the aforementioned Suess Effect, wherein d 13 C atm has changed in response to increased inputs of more isotopically negative CO 2 into the atmosphere (Keeling, 1979). The composition of CO 2 involved in the making of organic tissue is likely to affect the composition of that organic tissue (Arens, Jahren & Amundson, 2000); our study will test if this effect is compounded or mitigated by changes in other climate variables over this chronologically robust natural experiment.

Latitude
Latitude is expected to affect stomatal traits and therefore be related to carbon isotope fractionation due to its inverse relationship with light (specifically length of growing season and length of day) and temperature, and consequent effects on the maximum operating times for photosynthesis. A meta-analysis across 760 species in nine Chinese forest ecosystems showed a latitudinal variation in stomatal density and stomatal length at the community level (Wang et al., 2015). Given these morphological changes due to latitude, and the relationship between latitude and temperature, we might expect that D leaf values would be inversely related to latitude as well (Farquhar, Ehleringer & Hubick, 1989;Eq. (1)). The relationship between D leaf and latitude (15.9 S through 69.5 N) was observed in results of a meta-analysis of plants that used C 3 photosynthetic pathways (n ¼ 506) (Diefendorf et al., 2010), but any changes in D leaf as a function of latitude disappeared when latitude was decoupled from temperature and precipitation.

Precipitation
The stomata act as an inlet for CO 2 uptake as well as an outlet for leaf water loss via transpiration, which is why one might expect a relationship between carbon isotope fractionation and available water (represented in the paleo-record as reconstructed MAP). With decreased available water (i.e., decreased MAP) comes increased need for plant "water use efficiency" (as measured by the ratio of water used in photosynthesis to water lost through transpiration); plants therefore minimize water loss through the same stomata by fully closing, resulting in decreased carbon isotope fractionation (Farquhar, Ehleringer & Hubick, 1989). Previous meta-analyses, such as those by Diefendorf et al. (2010) and Kohn (2010), compared D leaf values of a wide variety of modern C 3 plants from many regions with MAP and found that D leaf varied significantly with MAP (p-value ¼ 0.0001 and R 2 ¼ 0.57; Diefendorf et al., 2010).

Temperature and seasonality
In addition to the climate variables with pre-established and applied relationships with D leaf values ([CO 2 ], MAP), various hypotheses have been proposed about the relationships between carbon isotope discrimination and other climate variables. For example, MAT could also constrain photosynthetic processes and associated carbon isotope fractionation because it gives a rough representation of extreme conditions and growing season length during which carbon assimilation occurs. For this reason, MAT is a well-addressed climate variable in previous isotope fractionation studies, but none have identified a relationship between fractionation and MAT (Arens, Jahren & Amundson, 2000;Diefendorf et al., 2010;Kohn, 2010;Schubert & Jahren, 2012). Furthermore, Helliker & Richter (2008) found that leaves maintained a constant internal temperature ideal for photosynthesis of 21.4 ± 2.2 C (total range of measurements), independent of external temperatures. In addition to MAT, maximum summer temperatures (particularly, the extreme highs associated with a warming climate) are expected to increase (Mirza, 2003). Increased maximum summer temperatures lead to increased evapotranspiration and more plant stress, which might affect carbon assimilation rates and stomatal conductivity (Farquhar, Ehleringer & Hubick, 1989;Diefendorf et al., 2010).
Seasonal variation is also thought to affect isotopic discrimination, resulting in a change in d 13 C leaf of up to 1-2‰ (Ehleringer, Phillips & Comstock, 1992;Arens, Jahren & Amundson, 2000), or in some deciduous trees such as maples, up to 6‰ between early spring and late fall (Lowdon & Dyck, 1974). Typically, more positive values are found in the winter (indicative of less isotopic discrimination) and more negative values occur in the summer (indicative of more discrimination). This effect strongest in arid and semiarid environments because they experience amplified seasonal temperature, precipitation, and evaporation effects (up to 4‰;Ehleringer, Phillips & Comstock, 1992). Our choice of sample locations should minimize this effect, as our specimens come from humid regions.
Though we do not expect a correlation between temperature and seasonality with D leaf values, nor are there good proxies for these variables in the fossil record, we include them here for completeness. This ensures that any observed noise is random or unaccountable for in the fossil record, rather than related to variables that are commonly measured and potentially relevant in plant isotope discrimination (Farquhar, Ehleringer & Hubick, 1989;Arens, Jahren & Amundson, 2000).

Finding a focal species
In addition to potentially confounding climate variables, variation in d 13 C leaf and D leaf values can be related to species-inherent carbon isotope fractionation. In a recent meta-analysis of C 3 plants conducted by Diefendorf et al. (2010) , 1992). However, the focus on geographic and climatic variability within that dataset resulted in a small number of analyses of any individual species. Thus, while some previous studies have proposed a universal D leaf value that represents C 3 plants on average (Arens, Jahren & Amundson, 2000;Gröcke, 2002), we focus here instead on an individual species (Thuja occidentalis; Cupressaceae) in order to avoid interspecific variation and phylogenetic/evolutionary effects in plant biochemistry.
Thuja occidentalis is a widespread evergreen gymnosperm with a distribution today extending throughout temperate deciduous and boreal forests in North America, and an extensive fossil record in localities across North America dating back to the Late Cretaceous (∼71 million years ago; LePage, 2003;Eckenwalder, 2009). T. occidentalis leaves have longer life spans (>1 year) than deciduous trees (Givnish, 2002;Pease, 1917), which makes them less vulnerable to seasonal variability and hardier in sedimentary archives (Diefendorf, Leslie & Wing, 2015). While some studies of other individual species have demonstrated unexplained internal isotopic variation of up to 3‰ (Tieszen, 1991), Mervenne (2015) found that D leaf values of T. occidentalis showed the least amount of internal isotopic variation within a single species grown in a common garden site (e.g., 18.91 ± 0.46‰ vs. Taxus: 20.05 ± 1.93‰) when compared to 56 species native to temperate forests. This makes T. occidentalis an excellent focal taxon for a single-species study.
This study incorporates the natural shifts in [CO 2 ] concentrations as driven by fossil fuel combustion and other anthropogenic inputs since the Industrial Revolution (280 ppm: Pre-Industrial, to ∼410 ppm in 2014; IPCC, 2014) and d 13 C atm values (from -6.5‰: Pre-Industrial to -8.5‰: present; Araus & Buxó, 1993;Elsig et al., 2009;White, Vaughn & Michel, 2015) to examine the relationship between T. occidentalis' carbon isotope fractionation and leaf chemistry (C:N ratios) within a range of climate variables. The patterns of D leaf values over a range of d 13 C atm values highlight the limitations of d 13 C leaf change as a tool for better understanding the biosphere and atmosphere. If d 13 C leaf values change in sync with d 13 C atm values, this could mean that carbon assimilation is continuing in this species as it was prior to Industrialization, perhaps indicating T. occidentalis' lack of adaptation to increased [CO 2 ]. Additionally, it may indicate that d 13 C leaf provides another way to track anthropogenic changes to the environment in the recent past and the future.

MATERIALS AND METHODS
We measured d 13 C leaf values of T. occidentalis extending from present-day to Pre-Industrial historical records using both newly collected and herbarium material. This included collecting leaf material of T. occidentalis specimens (n ¼ 142 collected between 1804 and 2017 with four unknown dates of collection) from across the Great Lakes region from herbaria (Figs. 1A and 1B; Table S1) and in natural present-day occurrences across a range of climate conditions (see Table S2). Thuja has small, 1-10 mm long scale-leaves adpressed along a small branch. Three cm portions of branches representing a single growth year with multiple scale-leaves were cleaned in an ultrasonic bath of deionized water to remove surface debris, oven dried at 50 C for 48 h and homogenized; this removes any within-leaf isotopic variation. Aliquots of each T. occidentalis specimen (0.6-0.8 mg) were placed into tin capsules and placed in a Costech elemental analyzer to measure %C and %N (as well as C:N ratio). Second aliquots of each T. occidentalis specimen (0.6-0.8 mg) were placed into tin capsules and placed in a combustion module inlet coupled to a Picarro G2201-i cavity ring-down spectrometer (CRDS) to measure d 13 C values of each specimen. Duplicates were run on both machines to insure homogeneity. Results of each CRDS run were internally calibrated using nine acetanilide standards (d 13 C ¼ -28.17 ± 0.16‰), two IAEA-600 caffeine standards (d 13 C ¼ -27.77 ± 0.04‰) and two IAEA-CH-6 sucrose To test whether other climate variables confounded the relationship between d 13 C atm and d 13 C leaf , we calculated carbon isotope fractionation values (D leaf ) for each specimen using known yearly d 13 C atm values (Eq. (1)). We derived global d 13 C atm values from direct (Rubino et al., 2013) and interpolated ice and firn core measurements according to the year of sampling (White, Vaughn & Michel, 2015). This does not account for ecosystem or microbiome-level deviations in d 13 C atm , but does account for the greater source isotopic value. D leaf values made it possible to de-couple the human-driven relationship between increased [CO 2 ] and decreased d 13 C leaf values (the Suess Effect; Keeling, 1979) because they were calculated by isolating leaf isotope ratios from the atmospheric isotopic signal.
For each specimen locality, environmental data expected to affect d 13 C leaf values (MAP, MAT, and maximum summer temperature; Farquhar, Ehleringer & Hubick, 1989) were derived from global databases using exact latitude and longitude coordinates of specimen origin (PRISM Climate Group, 2004;Fick & Hijmans, 2017;Government of Canada, 2018). All contiguous United States data was compiled from Oregon State University's PRISM database, which interpolates data from local weather stations at a resolution of four km (PRISM Climate Group, 2004). d 13 C leaf values were also compared with d 13 C atm and [CO 2 ] values at collection times, as retrieved from NOAA databases documenting values found at Mauna Loa Observatory in Hawaii (White, Vaughn & Michel, 2015) and measured on an isotope-ratio mass spectrometer at the institute of arctic and alpine research (INSTAAR) in the University of Colorado, Boulder. d 13 C leaf values for which d 13 C atm values were unavailable were not included in D leaf calculations or comparisons to climate variables (Table S3). While D leaf values combine d 13 C values measured in this experiment on the CRDS with d 13 C atm values measured on the IRMS at INSTAAR, we did not have access to individual errors for d 13 C atm and could not propagate the error. Therefore, we used our reproducibility error of 0.3‰, which is larger than the expected error for the IRMS (0.1‰), to be conservative.
Five-point moving averages of d 13 C leaf values were calculated to eliminate random noise caused by estimating older specimens' exact collection dates (a result of long collecting expeditions and limited recording resources). We regressed isotope values against climate variables (MAT, MAP, d 13 C atm and annual [CO 2 ]) to examine potential drivers of the d 13 C leaf values. We calculated the best-fit line using linear least squares regression to minimize the average distance between modeled y-values and actual y-values (d 13 C leaf and D leaf ) and calculated coefficients of determination. We used R 2 to determine predictive relationship between the given x-variable and y-variable. Additionally, we calculated p-values using the F-test to determine the chance of null hypothesis (p-value > 0.05).
The map of sampling location was created using R version 3.5.0 (R Core Team, 2014), and the ggplot2 (Wickham, 2016) and maps (v3.3.0, Becker & Wilks, 2018) packages. The full code is available in the Supplemental Files.

RESULTS
Values of d 13 C leaf ranged from -21.92‰ (collected in 1899) to -28.51‰ (collected in 2017), with a mean of -25.05 ± 1.32‰ (standard deviation; Table S2). D leaf values ranged from 15.11‰ to 20.97‰ (mean: 17.93 ± 1.11‰ standard deviation). Minimum, maximum, and mean values for climate variables are shown in Table S2. All data can be found in Table S3. D leaf vs. MAT, maximum summer temperature, latitude, seasonality, and MAP There was no relationship between D leaf values of T. occidentalis and MAT ( Fig. S1A; , nor between D leaf of T. occidentalis and maximum summer temperature ( Fig. S1B; R 2 ¼ 0.0051, p-value ¼ 0.51) in the temperature range listed in Table S2. No relationship was found between D leaf and latitude ( Fig. S1C; R 2 ¼ 0.0057, p-value ¼ 0.42). Additionally, no relationship was found using a multivariate linear regression approach to combine codependent variables: MAT and latitude (R 2 ¼ 0.002, p-values ¼ 0.92, 0.79 respectively). There was no relationship between D leaf of T. occidentalis and MAP ( Fig. 2; Fig. S1D; R 2 ¼ 0.0138, p-value ¼ 0.21), nor was there a relationship between D leaf of T. occidentalis and elevation ( Fig. S4A; R 2 ¼ 0.0138, p-value ¼ 0.55).

The atmosphere: d 13 C leaf and D leaf vs. [CO 2 ]
There was a linear relationship between d 13 C leaf and [CO 2 ] ( Fig. S3; R 2 ¼ 0.61, p-value < 0.001), but there was no relationship between D leaf and [CO 2 ] ( Fig. 4; Fig. S4; R 2 ¼ 0.0059, p-value ¼ 0.38). Because D leaf stayed constant, with simultaneous changes in d 13 C leaf and d 13 C atm , we can determine that changing d 13 C leaf with increased CO 2 was an effect of the changing isotopic composition of atmospheric CO 2 (d 13 C atm ) and was not related to [CO 2 ].   -4 robur, Fagus sylvatica, Cedrela odorata) where there was an influence. Additionally, regional variations in d 13 C atm resulting from proximity to respiring soil could provide a different baseline d 13 C atm value than the global one used based upon Mauna Loa Observatory's gas samples (Wehr & Saleska, 2015). While we were unable to control exactly where on the tree each sample was taken from nor the ecosystem-specific parameters that might influence d 13 C atm , each herbarium record includes notes on the sampling location, often including approximate maturity of the tree and height sampled from. Additionally, as expeditions are done without heavy machinery, it is likely that our specimens were sampled from approximately the height of a human, which is well out of the range of isotopic influence from soil respired CO 2 (Bazzaz & Williams, 1991). Finally, differences in d 13 C leaf values differ depending on where on the leaf isotopes are sampled from (Gao et al., 2015). We controlled for this effect by homogenizing several entire leaves per sample and running duplicate isotope analyses, which all came within machine error (0.3‰) of one another. Complex nitrogen dynamics, as well as those of other macro-and micro-nutrients, may play a role in vital effects that relate ultimately to carbon assimilation and/or carbon isotope fractionation, but these changes are very difficult and complicated to reconstruct in the geologic record and thus were not investigated in this study (Godfrey & Glass, 2011). It is not possible to understand shifts in regional nitrogen availability fully even within the historical record: without soil cores collected and preserved from the same sites and times of leaf collection, it is not straightforward to consider anything other than nitrogen content. Given that no one collected soils at each historical site, we cannot address changes in nitrogen dynamics quantitatively. Furthermore, for geologic applications fossil leaves and corresponding paleosols (fossil soils) are not typically preserved together. In order to account for the potential range in nitrogen dynamics, we focus here on using many specimens from across a wide landscape of different land-uses and thus a wide range of potential nitrogen dynamics, but under known [CO 2 ] conditions. To account for historical and geological limitations of reconstructing the nitrogen cycle, this study focuses instead on the relationship between systemic changes in C:N ratios in relation to d 13 C leaf , d 13 C atm , or as indication of changes in carbon dynamics.

D leaf , temperature, latitude and seasonality
There were a number of climate variables that we did not expect to have a correlative relationship with D leaf , but we addressed to ensure that they were not confounding variables (MAT, maximum summer temperature, latitude, seasonality). Due to the consistent internal temperature and lack of relationship between D leaf to MAT in previous studies, we expected no relationship between D leaf and MAT nor maximum summer temperature. As we expected, MAT showed no significant relationship with D leaf values neither in meta-analyses, growth chamber experiments (Diefendorf et al., 2010;Schubert & Jahren, 2012) nor this study (Fig. S1A). Because MAT and latitude are inherently related, we expected D leaf values to relate to latitude in the same way as when compared to MAT. As predicted, there was no relationship between D leaf of T. occidentalis and latitude. In order to ensure that these two codependent variables were treated as such, we ran a multiple linear regression, which gave a similarly low coefficient of determination (R 2 ¼ 0.002) comparing actual D leaf values with D leaf values predicted using this regression. Thus, we are confident that both treated as independent and co-dependent variables, MAT and latitude do not play a role in variance of D leaf values.
Additionally, we took into account variation in time of year specimen was collected and found no significant relationship between season and D leaf (Figs. S5A and S5B). This is likely because T. occidentalis is not deciduous and does not shed its leaves annually, thus, once homogenized, the isotopic composition of the leaf is representative of average discrimination during the leaf's exposure.

D leaf and [CO 2 ]
In this study we found no relationship between D leaf and [CO 2 ] from 280 to 410 ppm ( Fig. 4; Fig. S4). It is possible that the relationship observed by Schubert & Jahren (2012) exists only for idealized controlled growth chamber conditions and not in natural environments (Lomax et al., 2019). Alternatively, gymnosperms may respond more slowly than angiosperms to increases in [CO 2 ] due to their longer average lifespans and lack of senescence (Brodribb, Pittermann & Coomes, 2012). The changes in carbon assimilation as represented by increased carbon fractionation under short-term, ideal growth chamber conditions cannot be used to predict biological response to rapid changes in [CO 2 ]. In other words, plants, especially slow-growing woody plants, may not successfully adapt to anthropogenic changes of the present and future (Jump & Peñuelas, 2005).

D leaf and mean annual precipitation
When D leaf values of T. occidentalis were compared with MAP, Diefendorf et al.'s (2010) previously established relationship did not hold (Fig. 2), especially in low precipitation regimes (<1,000 mm yr -1 ) where the change in Diefendorf et al.'s D leaf values was most sensitive to changes in MAP. One explanation for the lack of relationship between D leaf of T. occidentalis and MAP is that this relationship breaks down on the single-species, or even plant functional type, level. In the aforementioned meta-analyses, plant functional type, species, and region were not controlled. D leaf values may be inherent to specific biomes but may not be representative of a general trend of any given plant or plant type to MAP. It is possible that the relationships seen in the meta-analyses by Diefendorf et al. (2010) and Kohn (2010) instead represent an array of taxon-specific constant isotopic values that collectively show a meta-relationship. This experiment could be further explored by performing the same experiment in natural settings across different biomes and different plant functional types. Assuming that D leaf values are indicative of water use (Givnish, 1979;Farquhar, Ehleringer & Hubick, 1989), this lack of relationship may also mean that plants with specific water use efficiencies and representative D leaf values are generally located in areas where they are not living in conditions that are stressed for water given their evolutionary adaptations. Geologically, this could mean that the presence of a particular taxon in fossil localities could provide a quantitative estimate for range of MAP, which could allow more specificity of paleoclimate regimes based on macrofossils (Nearest Living Relative and/Coexistence Approach; Mosbrugger & Utescher, 1997;Mosbrugger, 2009). In terms of future climate, this is indication that the chemistry of C 3 plants may not respond to regional changes as previously thought. This is of particular concern because the velocity of climate change, especially for continued high emission rate scenarios, is substantially faster than trees will be able to adapt to (Loarie et al., 2009;Diffenbaugh & Field, 2013).
It is also possible that the predicted relationship between D leaf and MAP is present for plants that are more responsive to their environment and/or have less extensive roots (to access deeper water sources), and thus the signal seen by Diefendorf et al. (2010) is a result of incorporating sensitive plants. Indeed, other studies using rapidly growing, highly sensitive herbaceous angiosperms have found a relationship between water treatment and carbon isotope discrimination (Lomax et al., 2019). However, the fossil record is biased toward preserving less sensitive, often woody plants due to preservation potential as well as presence within the fossil record (Looy et al., 2014); therefore, the utility of relationships based on highly sensitive plants may be muted in the fossil record.
Another explanation is that MAP is not an appropriate metric for measuring plant-available water, and while D leaf is still a measure of water use efficiency and this value is dynamic over conditions, snowmelt volumes and/or soil water-as driven by soil porosity and other factors-are better indicators of plant-available water. Further investigations using these variables will better constrain which water-related variables affect D leaf values of leaf tissues. The weak R 2 value between D leaf of T. occidentalis and MAP (0.0138) means that D leaf of T. occidentalis cannot be used to reconstruct paleo-MAP using the relationship determined by Diefendorf et al. (2010). Additional single-species experiments, particularly within angiosperms, should be conducted to look for correlations between D leaf and MAP to test whether the lack of relationship is due to a difference inherent to gymnosperms.
Carbon biomass (%C) and elemental leaf chemistry (%N, C:N ratios) as related to climate variables In addition to the response of d 13 C leaf values to climate variables, %C alone has been shown to respond directly to elevated [CO 2 ]. C i /C a ratios (the ratio of internal [CO 2 ] to atmospheric [CO 2 ]) of old growth T. occidentalis trees along Lac Duparquet, Quebec, increased under anthropogenic CO 2 fertilization, indicating tree response to enhanced CO 2 (Giguère-Croteau et al., 2019). This increase in tree productivity was demonstrated in the results of free-air concentration enrichment (FACE) experiments as well; in northern USA mid-latitude forests with loblolly pines (P. taeda), FACE experiment results indicated that elevated CO 2 induced increased carbon assimilation, resulting in increased carbon biomass, in woody tissues and increased %C of foliar storage as compared to trees grown under ambient CO 2 (Oren et al., 2001;Ainsworth & Long, 2005;Talhelm et al., 2013). Preliminary work in herbarium leaves found that increased [CO 2 ] related to Industrialization resulted in an increase in foliar %C with no change in %N (as source of N remains constant), and thus increased C:N ratios in some species (Mervenne, 2015). In order to contextualize changes in d 13 C leaf , this study examined coeval trends in leaf chemistry through elemental analysis of C and N (using N as a comparison point to see whether %C changes significantly with time).
These FACE experiments also showed that when run over longer time scales, trees reached a point of CO 2 acclimation and stopped increased carbon assimilation under enhanced CO 2 ; thus, predicted shifts in tree C-uptake may be short-lived, a pattern that will be inevitably discernable in a long-term study incorporating pre-Industrial leaf tissues through the present (Nowak, Ellsworth & Smith, 2004). Based on FACE experiments, we expected %C to have increased in leaves sampled from the early 1800s to the present, though we might see the rate of increase slow with time. However, we saw no relationship between %C, nor C:N ratios and time nor increase in [CO 2 ].
In fact, T. occidentalis specimens collected between 1804 and 2017 did not show changes in assimilation rates due to elevated CO 2 . C:N values of T. occidentalis showed no response to changes over time (with increased [CO 2 ]) or with atmospheric isotopic value (Figs. S2A and S2B). Though other organs in previous experiments responded to [CO 2 ], leaves, which are instrumental in the photosynthetic process as they are the organs directly in-taking atmospheric CO 2 , do not. A better understanding of all plant organ behavior is imperative to defining and quantifying potential carbon sinks or plant chemistry responses to global change (Goodale et al., 2002). d 13 C leaf and d 13 C atm Strong relationships have been found between above ground tissue and d 13 C atm values (p < 0.001; Jahren, Arens & Harbenson, 2008;Fig. 5), and this study provides a higher resolution look at the relationship between d 13 C leaf and d 13 C atm in a long-lived species within a natural system. In this initial natural experiment, the d 13 C leaf of T. occidentalis tracked changes in d 13 C atm (R 2 ¼ 0.74, p-value < 0.0001), mostly unencumbered by other climate factors. The slope for the linear relationship between d 13 C leaf and d 13 C atm is close to, but not exactly, 1, likely because the rate of change for d 13 C atm has not been linear, and acceleration in the change of d 13 C atm may not have been recorded immediately. Additionally, while there is no statistically significant relationship between any of the climate variables we tested and D leaf , it is unlikely that climate variables, especially in aggregate, play no role in carbon isotope discrimination within this species. Because d 13 C leaf showed a strong coefficient of determination with d 13 C atm , and no climate variables showed significant relationships with D leaf values, we can assume that d 13 C leaf values of modern T. occidentalis are strongly affected by d 13 C atm values. Additional work must be done to evaluate error in paleo uses of d 13 C leaf values of T. occidentalis, and future experiments should recreate more geologically reasonable conditions and climate changes (independent of anthropogenic factors). The relationship between d 13 C leaf and d 13 C atm values has implications for paleoclimate reconstructions of d 13 C atm as well as reconstructions of [CO 2 ] (Cerling et al., 1991;Franks et al., 2014). We emphasize how important it is to identify the value of d 13 C atm , such as in Tipple, Meyers & Pagani (2010) study, rather than just using the Pre-Industrial value of -6.5‰ (Cerling et al., 1991) because the d 13 C atm value has such a dramatic effect on the terrestrial part of the carbon cycle.

CONCLUSIONS
Though d 13 C leaf and D leaf values have been proposed as a proxy for [CO 2 ] and MAP based on previous research, this natural-world, species-controlled study shows no indication of such relationships. Thus, the use of D leaf values to reconstruct MAP and [CO 2 ] in the fossil record without taxonomic identification should be reconsidered. The relationship between d 13 C leaf and d 13 C atm values is more informative, and may provide a new proxy (d 13 C leaf values of Thuja) for reconstructing paleo-d 13 C atm or may indicate a lag in plant adaptation to unprecedentedly rapid climate change. Thuja extends up to 100 million years back to the Late Cretaceous, which makes this relationship potentially useful throughout the Cenozoic and into the Mesozoic era (Berry, 1915).
While this study focuses on one single species, further work is needed to assess other taxa at the species, genus, and family levels to examine whether the relationship between d 13 C atm and d 13 C leaf is consistent, and furthermore, generalizable. d 13 C leaf values of individual fossil leaves (in particular of Thuja leaves) cannot be used to reconstruct paleo-MAP as proposed by Kohn (2010), but average d 13 C leaf values of sites, as recorded in bulk soil organic matter, may allow us to predict precipitation ranges. Aboveground d 13 C leaf is thought to translate directly into the isotopic value of soil carbon (d 13 C org ; Arens, Jahren & Amundson, 2000). Bulk soil organic matter (d 13 C org ) is the combination of d 13 C of all decaying material from the ecosystem, with leaves especially abundant due to sheer volume. The average d 13 C leaf value of all trees found in a certain region will be found in the soil; therefore, soil d 13 C org values could be more reflective of particular precipitation at time of deposition than d 13 C leaf values. Further studies could evaluate the reliability of d 13 C org as a tool for MAP prediction and reconstruction.
This study implies constant carbon and nitrogen use and isotope fractionation relative to d 13 C atm by T. occidentalis. Due to the unprecedentedly rapid changes d 13 C atm and [CO 2 ] throughout Industrialization, this lack of change in carbon assimilation patterns, despite previous studies using modern d 13 C leaf values to reconstruct [CO 2 ], may indicate that modern systems are not appropriate analogues for many periods of the geologic record during which climate evolved more slowly. Modern climate change may be too rapid for plants to adapt, though more research should be done to evaluate whether this response is replicable in other species, genera, and plant functional types. It is possible that the pace of anthropogenic climate change makes modern relationships inappropriate analogues for paleoclimate.

Grant Disclosures
The following grant information was disclosed by the authors: The Herman and Margaret Sokol Foundation Award (University of Michigan). The Scott Turner Graduate Research Grant (University of Michigan).