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 ([CO2]) and isotopic value (δn class="Chemical">13Catm) of atmospheric CO2 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 CO2 (δ13Catm) 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. δ13Catm values provide a useful way to see changes in CO2 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 CO2 sources (e.g., methane, volcanism). δ13Catm values are particularly useful because they are parameters in models that reconstruct past changes to atmospheric [CO2] using paleosol carbonates (Cerling et al., 1991, 1992) or atmospheric [CO2] using plant stomatal parameters (Franks et al., 2014). Direct measurements of δ13Catm values only go back 50 years due to technological limitations, and longer-reaching ice core CO2 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 δ13Catm values in their own leaf n class="Chemical">carbon isotope values (δ13Cleaf) and fractionation values (Δ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 δ13C 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 (intercellular [CO2]) to C (atmospheric [CO2]), a ratio that is often used to represent water use efficiency.

While a and b are thought to be constant, we know that δ13Catm and C are changing rapidly. This could result in a corresponding change in Δleaf values as plants adapt to increased [n class="Chemical">CO2] or subsequent regional climate changes, for example, systematic changes in local precipitation). Alternatively, Δleaf values of leaves may stay constant but show marked changes in δ13Cleaf values corresponding to changes in δ13Catm values if leaves are incorporating δ13Catm into leaf tissues at a rate unaffected by other climate conditions.

Carbon Isotopes Related to Climate Variables

Previous studies have related Δleaf values to climate variables such as MAP, n class="Chemical">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), [CO2] (Schubert & Jahren, 2012, 2018), altitude (Korner, Farquhar & Wong, 1991; O’Leary, 1993), seasonality (Ehleringer, Phillips & Comstock, 1992), and δ13Catm 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 ([CO2] and δ13Catm 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

[CO2] and elevation

We would expect higher [CO2] to affect biochemical discrimination because of the known effects elevated [n class="Chemical">CO2] 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 [CO2] with a significant (21%) decrease in stomatal conductance (the rate of passage of atmospheric CO2 into plant tissue; Medlyn et al., 2001). Increased [CO2] 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 [CO2] (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 Δleaf on [CO2] may in part be due to isotopic discrimination associated with photorespiration (Schubert & Jahren, 2018). Indeed, previous growth chamber studies in prescribed CO2 environments showed increased carbon isotope fractionation with increased [CO2]; Schubert & Jahren (2012) found a strong hyperbolic correlation (r > 0.94) between [CO2] and Δleaf values in two species of herbaceous angiosperms. Based upon these growth chamber experiments (Schubert & Jahren, 2012, 2018), the relationship between [CO2] and Δleaf is expected to be most sensitive at geologically low [CO2] (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 clin class="Disease">mate variables such as temperature, vapor pressure, partial pressure of CO2 (pCO2), soil [CO2], 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 δ13Cleaf with increased altitude without controlling for other climate variables (Van de Water, Leavitt & Betancourt, 2002). In 2010, a meta-analysis assessing carbon isotope fractionation and discrimination values across a wide range of C3 plants (Diefendorf et al., 2010) found that when combined with MAP, elevation explained 61% of variability in Δ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 [CO2] and δ13Catm.

δ13Catm

With the current increase in [CO2], we have observed the aforementioned Suess Effect, wherein δn class="Chemical">13Catm has changed in response to increased inputs of more isotopically negative CO2 into the atmosphere (Keeling, 1979). The composition of CO2 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 n class="Chemical">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 Δleaf values would be inversely related to latitude as well (Farquhar, Ehleringer & Hubick, 1989; Eq. (1)). The relationship between Δleaf and latitude (15.9°S through 69.5°N) was observed in results of a meta-analysis of plants that used C3 photosynthetic pathways (n = 506) (Diefendorf et al., 2010), but any changes in Δleaf as a function of latitude disappeared when latitude was decoupled from temperature and precipitation.

Precipitation

The stomata act as an inlet for n class="Chemical">CO2 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 Δleaf values of a wide variety of modern C3 plants from many regions with MAP and found that Δleaf varied significantly with MAP (p-value = 0.0001 and R2 = 0.57; Diefendorf et al., 2010).

Temperature and seasonality

In addition to the climate variables with pre-established and applied relationships with Δleaf values ([n class="Chemical">CO2], 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 δ13Cleaf of up to 1–2‰ (Ehleringer, Phillips & Comstock, 1992; Arens, Jahren & Amundson, 2000), or in some den class="Chemical">ciduous 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 Δ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 δn class="Chemical">13Cleaf and Δleaf values can be related to species-inherent carbon isotope fractionation. In a recent meta-analysis of C3 plants conducted by Diefendorf et al. (2010), Δleaf values ranged from 13.4‰ for Pinus edulis in Utah, USA (Van de Water, Leavitt & Betancourt, 2002) to 28.4‰ for Cryptocarya concinna in Guangdong Province, China (Ehleringer, Phillips & Comstock, 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 Δleaf value that represents C3 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 den class="Chemical">ciduous 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 Δ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 [CO2] 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 δn class="Chemical">13Catm 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 Δleaf values over a range of δ13Catm values highlight the limitations of δ13Cleaf change as a tool for better understanding the biosphere and atmosphere. If δ13Cleaf values change in sync with δ13Catm 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 [CO2]. Additionally, it may indicate that δ13Cleaf provides another way to track anthropogenic changes to the environment in the recent past and the future.

Materials and Methods

We measured δ13Cleaf values of n class="Species">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 δ13C 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 (δ13C = −28.17 ± 0.16‰), two IAEA-600 caffeine standards (δ13C = −27.77 ± 0.04‰) and two IAEA-CH-6 sucrose standards (δ13C = −10.45 ± 0.03‰) in each run, as seen in Cotton, Sheldon & Strömberg (2012) study. Reproducibility of replicate analyses was better than 0.3‰.

Figure 1

Map of locations of Thuja occidentalis specimens.

Locations of (A) collections of Thuja occidentalis from across the world. (B) Specimens collected in eastern North America (in the typical habitat of Thuja occidentalis).

Locations of (A) collections of Thuja occidentalis from across the world. (B) Spen class="Chemical">cimens collected in eastern North America (in the typical habitat of Thuja occidentalis).

To test whether other climate variables confounded the relationship between δn class="Chemical">13Catm and δ13Cleaf, we calculated carbon isotope fractionation values (Δleaf) for each specimen using known yearly δ13Catm values (Eq. (1)). We derived global δ13Catm 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 δ13Catm, but does account for the greater source isotopic value. Δleaf values made it possible to de-couple the human-driven relationship between increased [CO2] and decreased δ13Cleaf 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 δn class="Chemical">13Cleaf 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). δ13Cleaf values were also compared with δ13Catm and [CO2] 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. δ13Cleaf values for which δ13Catm values were unavailable were not included in Δleaf calculations or comparisons to climate variables (Table S3). While Δleaf values combine δ13C values measured in this experiment on the CRDS with δ13Catm values measured on the IRMS at INSTAAR, we did not have access to individual errors for δ13Catm 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 δ13Cleaf values were calculated to eliminate random noise caused by estin class="Disease">mating older specimens’ exact collection dates (a result of long collecting expeditions and limited recording resources). We regressed isotope values against climate variables (MAT, MAP, δ13Catm and annual [CO2]) to examine potential drivers of the δ13Cleaf 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 (δ13Cleaf and Δleaf) and calculated coefficients of determination. We used R2 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 δ13Cleaf 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). Δleaf values ranged from 15.11‰ to 20.97‰ (mean: 17.93 ± 1.11‰ standard deviation). Minimum, maximum, and mean values for clin class="Disease">mate variables are shown in Table S2. All data can be found in Table S3.

Δleaf vs. MAT, maximum summer temperature, latitude, seasonality, and MAP

There was no relationship between Δleaf values of T. occidentalis and n class="Disease">MAT (Fig. S1A; R2 = 0.0152, p-value = 0.19), nor between Δleaf of T. occidentalis and maximum summer temperature (Fig. S1B; R2 = 0.0051, p-value = 0.51) in the temperature range listed in Table S2. No relationship was found between Δleaf and latitude (Fig. S1C; R2 = 0.0057, p-value = 0.42). Additionally, no relationship was found using a multivariate linear regression approach to combine codependent variables: MAT and latitude (R2 = 0.002, p-values = 0.92, 0.79 respectively). There was no relationship between Δleaf of T. occidentalis and MAP (Fig. 2; Fig. S1D; R2 = 0.0138, p-value = 0.21), nor was there a relationship between Δleaf of T. occidentalis and elevation (Fig. S4A; R2 = 0.0138, p-value = 0.55).

Figure 2

Carbon isotope fractionation values (Δleaf) vs. mean annual precipitation for Thuja occidentalis.

Δleaf vs. mean annual precipitation for Thuja occidentalis (black filled diamonds; Δleaf = −0.0003(MAP) + 18.75, R2 < 0.01). Data are compared with Diefendorf et al.’s (2010) global study (gray circles; Δleaf = 5.54(logMAP) + 4.07, R2 = 0.55), showing C3 specimens growing from 147 to 3,700 mm yr−1. Error bars along the y-axis represent the ±0.3‰ replicate reproducibility of standards.

Δleaf vs. mean annual precipitation for n class="Species">Thuja occidentalis (black filled diamonds; Δleaf = −0.0003(MAP) + 18.75, R2 < 0.01). Data are compared with Diefendorf et al.’s (2010) global study (gray circles; Δleaf = 5.54(logMAP) + 4.07, R2 = 0.55), showing C3 specimens growing from 147 to 3,700 mm yr−1. Error bars along the y-axis represent the ±0.3‰ replicate reproducibility of standards.

There was no relationship between Δleaf of T. occidentalis and month of collection. T-test results showed that the mean residual Δleaf values for each season were not significantly different (Spring: 0.19, Summer: −0.13, Fall: −0.01, Winter: −0.20; Fig. S5A; Month-by-month R2 = 0.0057; Fig. S5B).

%C, %N and C:N ratios

Values of %C ranged from 30.56 weight % to 61.44 weight % (±7.90%), with a dataset average of 48.86 %C. %N values ranged from 0.55 weight % to 2.28 weight %, with a dataset average of 1.31% (±0.21%). C:N ratios ranged from 20.1:1 (due to high %N) to 86.2:1 with an average C:N value of 40.7:1 (±7.9 C:N). There was no relationship between C:N ratios and time, δ13Cleaf, Δleaf, or [n class="Chemical">CO2] (R2 = 0.0227, p-value = 0.25; R2 = 0.0297, p-value = 0.68; R2 = 0.0157, p-value = 0.94; and R2 = 0.0029, p-value = 0.18, respectively; Fig. 3; Figs. S2A–S2C).

Figure 3

C:N ratios of specimens vs. [CO2].

C:N ratios of specimens vs. [CO2] (ppm) from Pre-Industrial values of 280–410 ppm. Error bars are associated with the 0.9% replicate reproducibility of standards.

C:N ratios of specimens vs. [n class="Chemical">CO2] (ppm) from Pre-Industrial values of 280–410 ppm. Error bars are associated with the 0.9% replicate reproducibility of standards.

The atmosphere: δ13Cleaf and Δleaf vs. [CO2]

There was a linear relationship between δ13Cleaf and [n class="Chemical">CO2] (Fig. S3; R2 = 0.61, p-value < 0.001), but there was no relationship between Δleaf and [CO2] (Fig. 4; Fig. S4; R2 = 0.0059, p-value = 0.38). Because Δleaf stayed constant, with simultaneous changes in δ13Cleaf and δ13Catm, we can determine that changing δ13Cleaf with increased CO2 was an effect of the changing isotopic composition of atmospheric CO2 (δ13Catm) and was not related to [CO2].

Figure 4

Carbon isotope fractionation values (Δleaf) vs. [CO2] for Thuja and growth chamber dicots.

Δleaf vs. [CO2] for Thuja and dicots. Shown are Thuja occidentalis (black empty triangles; Δleaf = 0.0029([CO2]) + 17.04, R2 = 0.0085) and Schubert & Jahren’s (2012) growth chamber studies with replicates of Raphanus sativus and Arabidopsis thaliana (black and grey smooth lines). Schubert & Jahren’s data shows C3 specimens growing at 15 levels of [CO2] from 370 to 4,200 ppm.

Δleaf vs. [CO2] for Thuja and n class="Species">dicots. Shown are Thuja occidentalis (black empty triangles; Δleaf = 0.0029([CO2]) + 17.04, R2 = 0.0085) and Schubert & Jahren’s (2012) growth chamber studies with replicates of Raphanus sativus and Arabidopsis thaliana (black and grey smooth lines). Schubert & Jahren’s data shows C3 specimens growing at 15 levels of [CO2] from 370 to 4,200 ppm.

δ13Catm and δ13Cleaf

There was a relationship between δ13Cleaf and δn class="Chemical">13Catm (Fig. 5; R2 = 0.74, p-value < 0.001) as represented by Eq. (2). The y-intercept (−16.52) represents the average offset between δ13Catm and δ13Cleaf, otherwise expressed as the fractionation value, Δleaf. The slope of 1.20 (with 95% confidence intervals between 1.07 and 1.32) further indicates relatively little impact of other environmental variables on leaf fractionation from atmospheric CO2. This relationship was compared to that as extrapolated from the regression found in Jahren, Arens & Harbenson (2008; Eq. (3); Fig. 5) study using Raphanus sativus in growth chamber experiments.

Figure 5

δ13Cleaf vs. δ13Catm values of Thuja occidentalis.

There is a linear relationship between δ13Cleaf and δ13Catm values (‰) of T. occidentalis, defined by δ13Cleaf = 1.20 (±0.06 standard error) * (δ13Catm) − 16.52(±0.44 standard error) with an R2 value of 0.74 and a p-value <0.001. Error bars along the y-axis represent the ±0.3‰ replicate reproducibility of standards, error bars along the x-axis represent the ±0.1‰ replicate reproducibility associated with the isotope-ratio mass spectrometer at the Institute of Arctic and Alpine Research (INSTAAR) in Boulder, Colorado. Data points represent specimens collected for this study. The regression line below is derived from Jahren, Arens & Harbenson (2008), which used Raphanus sativus grown under elevated [CO2] within growth chambers to look at the relationship between δ13Cleaf and δ13Catm values.

There is a linear relationship between δ13Cleaf and δn class="Chemical">13Catm values (‰) of T. occidentalis, defined by δ13Cleaf = 1.20 (±0.06 standard error) * (δ13Catm) − 16.52(±0.44 standard error) with an R2 value of 0.74 and a p-value <0.001. Error bars along the y-axis represent the ±0.3‰ replicate reproducibility of standards, error bars along the x-axis represent the ±0.1‰ replicate reproducibility associated with the isotope-ratio mass spectrometer at the Institute of Arctic and Alpine Research (INSTAAR) in Boulder, Colorado. Data points represent specimens collected for this study. The regression line below is derived from Jahren, Arens & Harbenson (2008), which used Raphanus sativus grown under elevated [CO2] within growth chambers to look at the relationship between δ13Cleaf and δ13Catm values.

Discussion

Potentially complicating factors

Herbarium specimens are a useful way to look in high-resolution on this time scale, but were not necessarily consistently sampled between expeditions and years. Thus some factors such as n class="Disease">maturity and height of tree, which have been shown to relate to carbon isotope discrimination (Brienen et al., 2017), are not specifically accounted for here. However, Brienen et al. (2017) found that height was not a significant factor in δ13Cleaf values of the gymnosperm used (P. sylvestris), in contrast to broadleaf species (Quercus robur, Fagus sylvatica, Cedrela odorata) where there was an influence. Additionally, regional variations in δ13Catm resulting from proximity to respiring soil could provide a different baseline δ13Catm 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 δ13Catm, 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 CO2 (Bazzaz & Williams, 1991). Finally, differences in δ13Cleaf 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 ultin class="Disease">mately 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 [CO2] 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 δ13Cleaf, δ13Catm, or as indication of changes in carbon dynamics.

Δleaf, temperature, latitude and seasonality

There were a number of climate variables that we did not expect to have a correlative relationship with Δleaf, but we addressed to ensure that they were not confounding variables (n class="Disease">MAT, maximum summer temperature, latitude, seasonality). Due to the consistent internal temperature and lack of relationship between Δleaf to MAT in previous studies, we expected no relationship between Δleaf and MAT nor maximum summer temperature. As we expected, MAT showed no significant relationship with Δ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 Δleaf values to relate to latitude in the same way as when compared to n class="Disease">MAT. As predicted, there was no relationship between Δ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 (R2 = 0.002) comparing actual Δleaf values with Δ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 Δleaf values.

Additionally, we took into account variation in time of year specimen was collected and found no significant relationship between season and Δleaf (Figs. S5A and S5B). This is likely because n class="Species">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.

Δleaf and [CO2]

In this study we found no relationship between Δleaf and [CO2] 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 n class="Species">angiosperms to increases in [CO2] 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 [CO2]. 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).

Δleaf and mean annual precipitation

When Δleaf values of T. occidentalis were compared with MAP, Diefendorf et al.’s (2010) previously established relationship did not hold (Fig. 2), espen class="Chemical">cially in low precipitation regimes (<1,000 mm yr−1) where the change in Diefendorf et al.’s Δleaf values was most sensitive to changes in MAP. One explanation for the lack of relationship between Δ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. Δ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 Δ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 Δ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 C3 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 Δ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 n class="Species">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 Δleaf is still a measure of n class="Chemical">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 Δleaf values of leaf tissues. The weak R2 value between Δleaf of T. occidentalis and MAP (0.0138) means that Δ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 Δ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 δ13Cleaf values to clin class="Disease">mate variables, %C alone has been shown to respond directly to elevated [CO2]. C/C ratios (the ratio of internal [CO2] to atmospheric [CO2]) of old growth T. occidentalis trees along Lac Duparquet, Quebec, increased under anthropogenic CO2 fertilization, indicating tree response to enhanced CO2 (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 CO2 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 CO2 (Oren et al., 2001; Ainsworth & Long, 2005; Talhelm et al., 2013). Preliminary work in herbarium leaves found that increased [CO2] 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 δ13Cleaf, 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 CO2 acclin class="Disease">mation and stopped increased carbon assimilation under enhanced CO2; 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 [CO2].

In fact, T. occidentalis spen class="Chemical">cimens collected between 1804 and 2017 did not show changes in assimilation rates due to elevated CO2. C:N values of T. occidentalis showed no response to changes over time (with increased [CO2]) or with atmospheric isotopic value (Figs. S2A and S2B). Though other organs in previous experiments responded to [CO2], leaves, which are instrumental in the photosynthetic process as they are the organs directly in-taking atmospheric CO2, 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).

δ13Cleaf and δ13Catm

Strong relationships have been found between above ground tissue and δ13Catm values (p < 0.001; Jahren, Arens & Harbenson, 2008; Fig. 5), and this study provides a higher resolution look at the relationship between δn class="Chemical">13Cleaf and δ13Catm in a long-lived species within a natural system. In this initial natural experiment, the δ13Cleaf of T. occidentalis tracked changes in δ13Catm (R2 = 0.74, p-value < 0.0001), mostly unencumbered by other climate factors. The slope for the linear relationship between δ13Cleaf and δ13Catm is close to, but not exactly, 1, likely because the rate of change for δ13Catm has not been linear, and acceleration in the change of δ13Catm may not have been recorded immediately. Additionally, while there is no statistically significant relationship between any of the climate variables we tested and Δleaf, it is unlikely that climate variables, especially in aggregate, play no role in carbon isotope discrimination within this species. Because δ13Cleaf showed a strong coefficient of determination with δ13Catm, and no climate variables showed significant relationships with Δleaf values, we can assume that δ13Cleaf values of modern T. occidentalis are strongly affected by δ13Catm values. Additional work must be done to evaluate error in paleo uses of δ13Cleaf values of T. occidentalis, and future experiments should recreate more geologically reasonable conditions and climate changes (independent of anthropogenic factors). The relationship between δ13Cleaf and δ13Catm values has implications for paleoclimate reconstructions of δ13Catm as well as reconstructions of [CO2] (Cerling et al., 1991; Franks et al., 2014). We emphasize how important it is to identify the value of δ13Catm, 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 δ13Catm value has such a dramatic effect on the terrestrial part of the carbon cycle.

Conclusions

Though δ13Cleaf and Δleaf values have been proposed as a proxy for [n class="Chemical">CO2] and MAP based on previous research, this natural-world, species-controlled study shows no indication of such relationships. Thus, the use of Δleaf values to reconstruct MAP and [CO2] in the fossil record without taxonomic identification should be reconsidered. The relationship between δ13Cleaf and δ13Catm values is more informative, and may provide a new proxy (δ13Cleaf values of Thuja) for reconstructing paleo- δ13Catm 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 spen class="Chemical">cies, genus, and family levels to examine whether the relationship between δ13Catm and δ13Cleaf is consistent, and furthermore, generalizable. δ13Cleaf values of individual fossil leaves (in particular of Thuja leaves) cannot be used to reconstruct paleo-MAP as proposed by Kohn (2010), but average δ13Cleaf values of sites, as recorded in bulk soil organic matter, may allow us to predict precipitation ranges. Aboveground δ13Cleaf is thought to translate directly into the isotopic value of soil carbon (δ13Corg; Arens, Jahren & Amundson, 2000). Bulk soil organic matter (δ13Corg) is the combination of δ13C of all decaying material from the ecosystem, with leaves especially abundant due to sheer volume. The average δ13Cleaf value of all trees found in a certain region will be found in the soil; therefore, soil δ13Corg values could be more reflective of particular precipitation at time of deposition than δ13Cleaf values. Further studies could evaluate the reliability of δ13Corg as a tool for MAP prediction and reconstruction.

This study implies constant carbon and n class="Chemical">nitrogen use and isotope fractionation relative to δ13Catm by T. occidentalis. Due to the unprecedentedly rapid changes δ13Catm and [CO2] throughout Industrialization, this lack of change in carbon assimilation patterns, despite previous studies using modern δ13Cleaf values to reconstruct [CO2], 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.

Δ values vs. climate variables.

a) Δ values vs. mean annual temperature (−5 to 14°C) for collected specimens (R2 = 0.0001). b) Δ values vs. maximum summer (growing season) temperature (20 to 33°C) for collected spen class="Chemical">cimens (R2 = 0.0005). c) Δ values vs. latitude (43 to 63°N) for collected specimens (R2 = 0.0002). d) Δ values vs. mean annual precipitation (288 to 1724 mm yr−1) for collected specimens (R2 = 0.0005).

Click here for additional data file.

C:N vs. climate and temporal variables.

C:N ratios as determined by weight %C and weight %N measured on a Costech Elemental Analyzer compared to a) time collected (year), b) δ13Catm values (‰), c) δn class="Chemical">13Cleaf values (‰). Error bars on the y-axis are associated with the 0.9% replicate reproducibility of standards. Error bars on the x-axis in c) are associated with the 0.3‰ replicate reproducibility of the Picarro CRDS standards.

Click here for additional data file.

δ13Cleaf values vs. pCO2.

Linear regression between δ13Cleaf values (‰) and [n class="Chemical">CO2] (283 to 407 ppm) for collected specimens (R2 = 0.61). Error bars along the y-axis represent the ±0.3‰ replicate reproducibility of standards.

Click here for additional data file.

Δ vs. elevation & pCO2.

a) Linear regression between Δ values vs. elevation (4 to 1617 m above sea level). Error bars along the y-axis represent the ±0.3‰ replicate reproducibility of standards b) Δ and n class="Chemical">pCO2 (283 to 407 ppm) for collected specimens (R2 = 0.01). Error bars along the y-axis represent the ±0.3‰ replicate reproducibility of standards.

Click here for additional data file.

δ13Cleaf and seasonality.

a) Box and Whisker plot showing median (lines), mean and range of spring, summer, fall and winter residual values of δ13Cleaf. T-tests assuming unequal variance comparing Spring vs. Summer, Fall, Winter, Summer vs. Fall, Winter, and Fall vs. Winter could not prove a null hypothesis; that the mean for all seasonal collections was significantly different. Seasons were categorized by mean monthly temperatures <°2C, from 2–15°C and rising, >15°C, and between 2–15°C and falling in the Great Lakes Region (centralized lower Peninsula Michigan), b) Scatter plot including residual from mean δn class="Chemical">13Cleaf values as divided by month.

Click here for additional data file.

Click here for additional data file.

The range and mean values for measured climate variables in this study.

The minimum, maximum, and mean values for climate conditions for each spen class="Chemical">cimen as obtained from NOAA ESRL Global Monitoring Division (2016), White et al. 2015, and PRISM Climate Group, as well as Government of Canada (Ed.). (2018, January 11).

Click here for additional data file.

Raw Data.

*Direct Mauna Loa measurements (NOAA ESRL Global Monitoring Division (2016).

**White et al. (2015)) and Petit et al. (2001) ice core records.

***: for United States Climate Data: PRISM Clin class="Disease">mate Group Oregon State University. (2017). Data Explorer: Time Series Values for Individual Locations. Retrieved from Northwest Alliance for Computational Science & Engineering database. For Canadian Climate Data: Government of Canada (Ed.). (2018, January 11). All specimens used for isotope analysis are stored in the University of Michigan’s Earth Systems Laboratory (Dr. Nathan Sheldon & Dr. Ingrid Hendy) within the University of Michigan Earth and Environmental Science Department.

Click here for additional data file.

Click here for additional data file.