Tree-ring isotopes suggest atmospheric drying limits temperature-growth responses of treeline bristlecone pine.

Altitudinally separated bristlecone pine populations in the White Mountains (California, USA) exhibit differential climate-growth responses as temperature and tree-water relations change with altitude. These populations provide a natural experiment to explore the ecophysiological adaptations of this unique tree species to the twentieth century climate variability. We developed absolutely dated annual ring-width chronologies, and cellulose stable carbon and oxygen isotope chronologies from bristlecone pine growing at the treeline (~3500 m) and ~200 m below for the period AD 1710-2010. These chronologies were interpreted in terms of ecophysiological adaptations to climate variability with a dual-isotope model and a leaf gas exchange model. Ring widths show positive tree growth anomalies at treeline and consistent slower growth below treeline in relation to the twentieth century warming and associated atmospheric drying until the 1980s. Growth rates of both populations declined during and after the 1980s when growing-season temperature and atmospheric vapour pressure deficit continued to increase. Our model-based interpretations of the cellulose stable isotopes indicate that positive treeline growth anomalies prior to the 1980s were related to increased stomatal conductance and leaf-level transpiration and photosynthesis. Reduced growth since the 1980s occurred with a shift to more conservative leaf gas exchange in both the treeline and below-treeline populations, whereas leaf-level photosynthesis continued to increase in response to rising atmospheric CO2 concentrations. Our results suggest that warming-induced atmospheric drying confounds positive growth responses of apparent temperature-limited bristlecone pine populations at treeline. In addition, the observed ecophysiological responses of attitudinally separated bristlecone pine populations illustrate the sensitivity of conifers to climate change.

Introduction

Tree-ring archives from the unique and ancient bristlecone pine (Pinus longaeva D. K. Bailey) contain annually resolved information on historic changes in temperature and precipitation covering several millennia (Ferguson 1968, LaMarche and Stockton 1974, Hughes and Funkhouser 1998, Salzer, Bunn, et al. 2014). A remarkable aspect of this species’ climate–growth response is that an altitudinal separation of as little as 80 vertical metres can differentiate nearby populations in terms of apparent temperature limitation (near treeline) and moisture limitation (at lower elevations) (Salzer et al. 2009, Tran et al. 2017). Tree-ring chronologies from bristlecone pine growing in the White Mountains of California (USA) revealed that trees growing at the modern treeline (situated at ~3500 m above mean sea level) exhibit positive growth anomalies during most of the twentieth century, while trees growing at lower elevations show little changes in growth rate over this time period (Lamarche 1974, Salzer et al. 2009). Treeline populations on south-facing slopes have shown a reversal of the earlier positive growth anomalies since the mid-1990s, which suggests the onset of moisture limitation at higher elevations (Salzer, Larson, et al. 2014). However, the specific ecophysiological mechanisms underlying these twentieth century growth responses remain to be determined (Tran et al. 2017).

Growth at the treeline is generally considered to be limited by temperatures suboptimal for xylogenesis (i.e., the formation of water conducting vascular tissue) (Rossi et al. 2007, Ziaco, Biondi, Rossi, et al. 2016), although site-specific and species-specific factors often influence growth as well (Körner and Paulsen 2004, Holtmeier and Broll 2005, Barbeito et al. 2012, Greenwood et al. 2015). Moisture limitation becomes more prevalent at lower elevations owing to higher temperatures and increased atmospheric vapour pressure deficit (VPD) (Grace et al. 2002, Wiley and Helliker 2012, Körner 2015). Mechanistically, water limitation reflects drought stress as a result of hydraulic constraints on plant gas exchange (Sperry 2000). Drought stress specifically stems from too low water potentials in the soil-leaf continuum as a consequence of limited soil water supply and/or high atmospheric demand (Venturas et al. 2017). Too low xylem water potential can inhibit xylogenesis directly by offsetting cellular metabolic activity and cell turgor (Lautner 2013, Cuny and Rathgeber 2016) and indirectly as drought-induced stomatal closure limits photosynthesis and constrains the supply of carbohydrates required for growth (Sevanto et al. 2014, Adams et al. 2017, Choat et al. 2018).

Recent work highlights the role of tree hydraulics in modulating the climate–growth responses in coniferous species (Castagneri et al. 2015, McDowell et al. 2016, Prendin et al. 2018). Conifers are particularly sensitive to changes in VPD because this group maintains relatively wide hydraulic safety margins in order to minimize embolisms in the stem (Johnson et al. 2012). As a result, conifers generally reveal earlier drought-induced stomatal closure than angiosperms (Meinzer et al. 2009, Carnicer et al. 2013). Although this conservative evaporation strategy protects conifers from hydraulic failure in the stem, it comes at the expense of reduced leaf-level photosynthesis during drought (Timofeeva et al. 2017). As warming is typically associated with rising VPD (Cook et al. 2014), warming can increase drought stress, especially in water limited environments (Williams et al. 2013). Conversely, warming may improve tree-water status when xylogenesis is limited by suboptimal temperatures and sufficient water is available. Earlier evidence further suggests that rising atmospheric CO2 concentrations (ca) can stimulate growth in bristlecone pine (Lamarche et al. 1984, Graybill and Idso 1993). This carbon fertilization effect is most likely to occur under moisture-limited conditions as a result of water savings from CO2-induced stomatal closure (Fatichi et al. 2016). Still, the question if and how rising ca stimulates tree growth remains a topic of active research and debate (Körner 2015, Sleen et al. 2015, Girardin et al. 2016, Terrer et al. 2016, Guillemot et al. 2017).

Stable carbon isotopes (δ13C) and stable oxygen isotopes (δ18O) measured on tree ring cellulose can provide further insight into the changes in stomatal conductance, transpiration and photosynthesis that occur in response to variable growth conditions (Scheidegger et al. 2000, McCarroll and Loader 2004, Waterhouse et al. 2004, Grams et al. 2007, Guerrieri et al. 2017). The environmental factors and plant-physiological responses that influence cellulose δ13C and δ18O are generally well understood; however, especially the interpretation of cellulose δ18O in terms of plant-physiological responses is difficult as changes in source water and atmospheric δ18O influence cellulose δ18O as well (Roden and Siegwolf 2012, Barbour and Song 2014). In theory, cellulose δ18O can reflect changes in VPD, stomatal conductance and transpiration owing to their combined influence on oxygen fractionation in transpiring leaves (Barbour et al. 2004, Barbour 2007, Farquhar et al. 2007). However, changes in atmospheric circulation patterns, soil evaporation processes, seasonal variations in weather conditions and back-diffusion of atmospheric water vapour can confound the ecophysiological interpretation of δ18O signals recorded in tree ring cellulose (Zimmermann et al. 1967, Gat 1980, Roden and Siegwolf 2012, Treydte et al. 2014, Lehmann et al. 2018). In the case of bristlecone pine growing in the White Mountains, the δ18O signal recorded in tree ring cellulose predominantly reflects changes in leaf evaporative enrichment during the growing season (Berkelhammer and Stott 2009). The resulting relatively weak influence of temporal changes in source water composition is expressed in the long-term stability of δ18O signal in bristlecone pine tree ring cellulose (Bale et al. 2010). Hence, combined records of δ13C and δ18O measured on bristlecone pine tree rings can potentially help to elucidate the responses of leaf gas exchange in relation to variations in tree ring growth by applying the theoretical principles of the conceptual dual-isotope model (Scheidegger et al. 2000). However, given the inherent uncertainty conveyed in cellulose δ18O, the dual-isotope model should be interpreted with caution (Roden and Siegwolf 2012).

We aim to explore how combined changes in temperature, precipitation, VPD and ca affect growth and gas exchange of bristlecone pine growing near treeline and at lower elevations. To achieve this, we developed absolutely dated annual ring width, and cellulose δ13C and δ18O chronologies from bristlecone pine populations growing at the modern treeline and approximately 200 m below that span the period from AD 1710 to 2010. The combination of annual ring widths and cellulose stable isotopes allows, in principle, to independently reconstruct changes in stem growth and leaf gas exchange, which are not necessarily synchronized and may respond to different environmental variables at different time scales. Our population selection includes the proposed altitudinal shift from moisture sensitivity to temperature sensitivity in a relatively small spatial area (Salzer, Larson, et al. 2014, Tran et al. 2017). Data were analysed in a three-staged approach. First, the chronologies of tree ring widths were analysed for their association with key climate variables. Second, coeval changes in cellulose δ13C and δ18O were interpreted in terms of changes in leaf gas exchange following the theoretical principles of the dual-isotope model (Scheidegger et al. 2000, Grams et al. 2007). Third, we modelled the contribution of temperature, VPD and ca to changes in leaf-level photosynthesis, transpiration and stomatal conductance using a semi-empirical leaf gas exchange model (de Boer et al. 2016). This combination of methods allows us to interpret differential growth responses of bristlecone pine in relation to the twentieth century climate variability and changes in leaf gas exchange.

Materials and methods

Site characteristics

We obtained cellulose core samples from 17 Pinus longaeva D. K. Bailey (bristlecone pine) individuals distributed over three growth localities on Sheep Mountain in the Inyo-White Range, Mono County, eastern California (USA) (see Fig. S1 available as Supplementary Data at Tree Physiology Online). The summit of Sheep Mountain has an elevation of 3809 m (all elevations are expressed relative to mean sea level) with the treeline situated at ~3500 m. The geology of Sheep Mountain is moderately diverse, consisting of complexly folded and faulted volcanic rocks and metamorphosed sandstone, shale, limestone and dolomite (Hall, 1991). Within this diverse geology, the bristlecone pine stands predominantly occur on dolomite substrate patches owing to the favourable microclimate (lower soil temperatures and transpiration) resulting from the greater reflectivity (Fritts 1976). Trees in sites 1 A (five individuals sampled) and 1B (six individuals sampled) are situated at elevations between 3482 m and 3523 m at the modern treeline, trees in site 2 (six individuals sampled) grow at slightly lower elevations between 3293 m and 3338 m (see Table S1 available as Supplementary Data at Tree Physiology Online). The approximate 200 metre elevation difference between sites is sufficient to include a shift from moisture limitation to temperature limitation (Salzer et al. 2009, Salzer, Larson, et al. 2014, Tran et al. 2017). All sites have south-facing aspects with a slope of ~20° with trees growing on comparable substrate.

Sampling of tree ring cellulose and foliar material

Trees in site 1 A were sampled on the 5 August 2010, the trees in sites 1B and 2 were sampled between the 23 and 25 August 2011. Trees were selected following the conditions of USFS permit WMD110026T obliging that no trees with previous identification tags and/or less than 25 cm strip-bark were sampled. All trees were cored using a 12 mm internal diameter Haglöf increment borer following the methods described by Stokes and Smiley (1996). Cores were taken at breast height, between 1.3 m and 1.4 m above the ground (Stokes and Smiley 1996, Wong and Lertzman 2001). Each tree was cored on the side slope wherever possible to prevent sampling of distorted growth ring patterns. All trees were tagged with an engraved aluminium metal tag and given a unique ID with the prefix ‘SWAN’. Strip-bark trees were tagged on dead-wood. After sampling, each core was placed in a perforated PVC tube for transport. Circumference of each tree was measured for the Diameter at Breast Height (DBH). Small branches were sampled from four treeline and four below-treeline individuals for foliar analyses (Table S1).

Tree-ring chronologies

Cores were mounted, without glue, in custom-made wooden blocks and sanded with progressively finer grades of abrasive paper and measured with a precision of 0.01 mm using a Velmex measuring stage. TSAP-Win v4 (Time Series Analysis and Presentation program for Windows) was used to construct and cross-date ring-width chronologies using standard dendrochronological techniques (Stokes and Smiley 1996). Chronology validation was achieved using COFECHA (Holmes 1983, Grissino-Mayer 2001). The ring-width series were absolutely dated against independent chronologies (Graybill and Idso 1993, Salzer et al. 2009). Tree ages were estimated using the pith-estimation technique (Duncan 1989). The earliest dated tree ring as well as the estimate age of each cored tree are shown in Table S2, available as Supplementary Data at Tree Physiology Online. We note that ring-width chronologies are often expressed in terms of basal area increments (BAIs) to more accurately represent changes in tree growth (Biondi and Qeadan 2008). However, this measure is somewhat problematic to use for our samples as some of the sampled individuals had non-circular cross-sections owing to their strip-bark morphology (Table S1). Our interpretation of differential growth responses is therefore based on ring widths.

Carbon and oxygen stable isotopes

Annual isotope measurements were performed on pooled samples of α-cellulose from all cored trees from the treeline and below-treeline sites, respectively. To control for potential bias in the pooled samples we also performed decadal isotope measurements on non-pooled samples from five individuals of the treeline sites 1 A and 1B, and five individuals of the below-treeline site 2, respectively (see Fig. S2 available as Supplementary Data at Tree Physiology Online). Each core was sliced with a sharp scalpel to form thin wood slivers (~40 μm). As the rings were small with a relatively minor proportion of latewood, no attempt was made to separate earlywood and latewood. Standard methods were employed to isolate α-cellulose with the oxidation of lignin using acidified sodium chlorite and the subsequent hydrolysis of hemicelluloses (Loader et al. 1997, Rinne et al. 2005). The resulting α-cellulose was homogenized using a Hielscher ultrasonic probe (Laumer et al. 2009) and freeze-dried in a Thermo Savant ModulyoD freeze drier at −45 °C for at least 48 h. For carbon and oxygen isotope analysis, α-cellulose (0.30–0.35 mg) was weighed into silver foil capsules and pyrolyzed over glassy carbon at 1090 °C in a PDZ Europa ANCA GSL elemental analyzer interfaced to a PDZ Europa 20–20 stable isotope ratio mass spectrometer (Young et al. 2011, Loader et al. 2015). Isotope ratios are expressed as per mille deviations relative to VPDB (carbon) and VSMOW (oxygen) standards. For an internal Sigma-Aldrich α-cellulose standard, analytical precision was typically ±0.1‰ for carbon and ±0.3‰ for oxygen.

The records of tree ring cellulose were corrected to a pre-industrial standard atmospheric δ13C value of −6.4‰ (δ13Catm) for plotting following (McCarroll et al. 2009) for plotting:

Here, δ13Ccor and δ13Ccell represent the corrected and uncorrected tree ring cellulose values, respectively.

Climate data

The climate variables used to explore relationships with the ring-width chronologies and cellulose stable isotopes were obtained from the CRU-TS4.00 data product at a 0.5° spatial resolution (Harris et al. 2014) and the PRISM AN81m data product at a 4 km resolution (Daly et al. 2008, 2015). From the CRU-TS4.00 data product, we extracted the temperature averaged over the growing-season months July–September (T) and cumulative precipitation over the preceding September-to-August period (Pp). The growing-season VPD (VPDgs) was derived from the PRISM data product. These climate variables are most relevant for bristlecone pine growth in the White Mountains (Salzer et al. 2009). The CRU climate variables were extracted from the grid cell centred on the latitude/longitude 37.25°N and −118.25°E. This single grid cell covers all growth localities sampled in our study. From PRISM, we selected two grid cells centred on 37.5324°N and −118.1997°E, and 37.5151°N and −118.1934°E covering the treeline and below-treeline sites, respectively. For each grid cell, the minimum and the maximum monthly average VPD were extracted and normalized relative to the period 1981–2010. The resulting anomalies were averaged between the two grid cells to obtain a single monthly average VPD record. The months July–September where subsequently averaged to calculate the anomalies in VPDgs. All climate variables were extracted for the period from 1901 (start of the CRU data product) to 2010 (youngest tree ring).

We note that CRU data product was used to derive multi-decadal climate trends in temperature and precipitation as the CRU data product is specifically developed and validated to represent climate trends and anomalies (Harris et al. 2014). The PRISM data product is specifically developed to quantify the influence of topography onsite-specific (hydro)climate in high spatial revolution. The PRISM data product is, however, potentially sensitive to changes in meteorological station cover, which makes the VPD data less suitable to represent multi-decadal trends (see Daly et al. 2008, 2015). Moreover, instrumental measurements of atmospheric humidity are increasingly scarce before the 1960s (Elliott 1995). VPD variability prior to the 1960s should, therefore, be interpreted to primarily reflect changes in temperature.

Conceptual dual-isotope model

Coeval changes in cellulose δ13C and δ18O can, in theory, be interpreted in terms of plant-physiological responses in relation to climate variability following the conceptual dual-isotope model (Scheidegger et al. 2000, Barbour et al. 2004, Grams et al. 2007). However, the signals should be interpreted with caution considering the many inherent uncertainties, especially in terms of cellulose δ18O which may reflect changes in source water and atmospheric δ18O as well (Roden and Siegwolf 2012). Moreover, the original dual-isotope model does not explicitly consider the effect of rising ca on plant-physiological responses and subsequent changes in the ratio of intercellular to atmospheric concentrations of CO2 (ci/ca) that may occur independent of climate changes (Farquhar et al. 1989, Ehleringer and Cerling 1995, Franks et al. 2013). Hence, interpreting the twentieth century cellulose δ13C records in relation to climate variability requires considering the effects of climate-independent plant-physiological responses to rising ca. Plant physiological responses in terms of ci to rising ca are typically grouped in three distinct response strategies: (i) maintaining a constant c, (ii) maintaining a constant drawdown in CO2 (constant ca–ci) and (iii) maintaining a constant ci/ca (Saurer et al. 2004, McCarroll et al. 2009, Voelker et al. 2016). Of the three strategies, the constant ci/ca response reflects optimization of investments required for water transport and carbon assimilation as predicted by theory (Prentice et al. 2014), whereas the constant ca–ci response is biased towards lower than theoretical optimal water use efficiency and the constant ci response is biased towards higher than optimal water use efficiency. To facilitate interpretation of the cellulose δ13C-dervied ci/ca records in terms of either of these responses, we detrended the cellulose δ13C-dervied ci/ca records relative to constant ci/ca response strategy, as explained in the following.

First, we derived the time-integrated ratio of intercellular to atmospheric concentrations of CO2 (ci/ca) recorded for each growing season in the tree ring cellulose: where Δ13C is the carbon discrimination calculated from δ13Ccell and δ13Catm (Farquhar et al. 1989), and a (4.4‰) and b (28‰) are constants representing fractionation due to diffusion and carboxylation, respectively (Farquhar et al. 1982).

The time-integrated leaf-interior CO2 concentration () was calculated from the cellulose δ13C-derived ci/ca and the Mauna Loa annual-average atmospheric CO2 record (ca-ML) (Keeling and Whorf 2004):

Subsequently, we expressed the deviation in ci (termed ∆ci) from the baseline constant ci/ca response for the annual population-averaged records of cellulose δ13C as

Here, ∆ci, and ca-ML are variables that change with calendar year t. The (constant) parameters ci,PI and ca,PI represent the average δ13C-derived ci and Mauna Loa-based ca over the pre-industrial period 1750–1849, respectively. The third term in Eq. (4) thus represents the baseline constant ci/ca response of ci to rising ca expressed relative to the pre-industrial average ci/ca of each population. Hence, the ∆ci represents the observed deviation in ci from the (moving) baseline response at time t (see Fig. S3 available as Supplementary Data at Tree Physiology Online). We note that this approach to detrending the cellulose δ13C-derived ci/ca records does not imply an a priori assumption on how plants adjust ci/ca to rising ca in the absence of climate change. Rather, the detrended records facilitate the interpretation of δ13C in terms of plant-physiological responses to rising ca that reflect one of the three typical strategies; no change in ∆ci means a constant ci/ca response consistent with optimization theory (Prentice et al. 2014), declining ∆ci reveals a bias towards higher water use efficiency, and rising ∆ci reveals a bias towards lower water use efficiency relative to the theoretical optimum.

Our interpretation of stable carbon and oxygen isotopes recorded in tree ring cellulose is subsequently based on coeval changes in ∆ci and δ18O. Changes in ∆ci are associated with the same-signed changes in stomatal conductance to CO2 (∆gsc) and photosynthesis (∆A), and opposite-signed changes in photosynthetic capacity (∆Amax), which are all expressed relative to the baseline constant ci/ca-response (as indicated by the ∆). Changes in cellulose δ18O are interpreted in terms of same-signed changes in VPD and/or opposite-signed changes in stomatal conductance to water vapour (gsw) and leaf-level transpiration (Barbour et al. 2004, Farquhar et al. 2007). For example, for a reduction in stomatal conductance in response to rising VPD, we expect to observe a reduction in ∆ci and an increase in δ18O. Owing to the proportionality between stomatal conductances to water vapour and CO2 (gsw ≈ gsc · 1.6), the dual-isotope model can logically rule-out potential conflicting interpretations of changes in leaf gas exchange parameters.

Modelling of leaf gas exchange

Considering the inherent uncertainty associated with the dual-isotope model (Roden and Siegwolf 2012), we also used a semi-empirical leaf gas exchange model (de Boer et al. 2016) to explore the responses of photosynthesis, transpiration and stomatal conductance to changes in temperature and VPD independent of cellulose δ18O. In short, this model used the annual records of cellulose δ13C-derived leaf interior partial pressure of CO2 () in combination with instrumental record of temperature to calculate the time-integrated photosynthesis at the leaf level () with the Farquhar et al. (1980) model of photosynthesis. Stomatal conductance and the time-integrated transpiration rate () were calculated with Fick’s law of diffusion and the instrumental records of VPD and ca. Our application of this model assumed that photosynthesis is generally limited by the rate of carboxylation. Values of the maximum carboxylation capacity normalized at 25 °C (Vcmax25) were estimated from the emperical relationship between leaf nitrogen content (Narea) and Vcmax25, following Walker et al. (2014) (see Table S3 and Fig. S4 available as Supplementary Data at Tree Physiology Online). The adjustment of Vcmax in response to changes in ca was modelled using the optimality approach of Wang et al. (2017), which reproduced the commonly observed down-regulation of Vcmax25 under rising CO2 (Ainsworth and Rogers 2007). The combination of methods involving the Farquhar model of photosynthesis, Fick’s law of diffusion and estimation of Vcmax25 from leaf nitrogen content is well established (Hikosaka et al. 2016), although changes in photosynthetic nitrogen use efficiency may occur in response to rising CO2 (Nijs et al. 1995). Uncertainty remains in environmental boundary conditions, especially considering the temporal variability in VPD prior to the 1960s (Elliott 1995). Details on the modelling approach are presented in the Supplemental Methods, available as Supplementary Data at Tree Physiology Online. An analysis of the model behaviour and sensitivity is presented by de Boer et al. (2016).

Model boundary conditions were based on a combination of site-specific modern temperature and VPD climatology (1981–2010) from the PRISM Norm81m data product at an 800 m resolution (Daly et al. 2008, 2015), temperature trends derived from the CRU-TS4.00 data product at a 0.5° spatial resolution (Harris et al. 2014), and VPD trends derived from the PRISM AN81m data product at an 4 km (Daly et al. 2008, 2015). To develop the model boundary conditions, we first used the PRISM Norm81m data product to calculate site-specific climatology for the period 1981–2010 in terms of minimum and maximum monthly temperature and VPD for the growing-season months July–September. The grid cell centred on the latitude/longitude 37.5151°N and −118.1934°E was used for the below-treeline locality, and the grid cell centred on the latitude/longitude 37.5324°N and −118.1997°E was used for the treeline localities. The PRISM site-specific climatology represents a temperature lapse rate ranging between 0.62°C and 0.83°C per 100 m, depending on the time of year. These lapse rates are slightly higher than standard adiabatic lapse rates (Van de Ven et al. 2007), but slightly lower than measured on the nearby Sheep Mountain during the growing season (Salzer, Bunn, et al. 2014). Second, we used the CRU-TS4.00 data product to obtain monthly temperature anomalies and the PRISM AN81m data product to obtain monthly VPD anomalies for the period 1901–2010 for each growing-season month. Anomalies were calculated relative to the 1981–2010 average to match the PRISM Norm81m climatology. The CRU-TS4.00 temperature was extracted from the grid cell centred on the latitude/longitude 37.25°N and −118.25°E. The PRISM VPD anomalies were calculated as an average from the two grid cells centred on the latitude/longitude 37.5324°N and −118.1997°E, and 37.5151°N and −118.1934°E covering the treeline and below-treeline sites, respectively. The monthly temperature and VPD anomalies were added to the minimum and maximum monthly average temperature and VPD from the site-specific PRISM Norm81m climate normals. The model boundary conditions for the treeline and below-treeline sites thus reflect common trends in growing-season temperature and VPD for the period 1901–2010 based on the CRU-TS4.00 and PRISM AN81m data products as well as site-specific differences in absolute growing-season temperature and VPD based on the high-resolution PRISM data representative of the period 1981–2010 (see Fig. S5 available as Supplementary Data at Tree Physiology Online).

With these boundary conditions, we performed four different simulations (Table 1): (i) the control ‘CO2, T & VPD’ simulation with observed ca, CRU-derived trend in site-specific temperature and PRISM-derived trend in VPD as boundary conditions, (ii) the ‘CO2’ simulation with observed ca together with constant early-20th century (1901–1930) site-specific temperature and constant VPD, (iii) the ‘CO2 & T’ simulation with observed ca, the CRU-derived trend in site-specific temperature, and constant early-twentieth century VPD and (iv) the ‘CO2 & VPD’ simulation used observed ca, the CRU-derived trend in site-specific VPD, and constant early-twentieth century temperature. The effect of rising ca was included in all simulations because the semi-empirical model requires input in terms of δ13C-derived ci/ca. This input variable cannot meaningfully be considered independent of changes in ca. The influence of temperature on photosynthesis was included via the temperature dependency of Rubisco-limited photosynthesis (Bernacchi et al. 2001).

Table 1.

Overview of the climate boundary conditions used in the leaf gas exchange modelling covering the period 1901–2010.

	Model boundary conditions
Simulation name	CO2	Temperature	VPD
CO2, T & VPD	Observed ca	CRU anomalies with PRISM climatology	PRISM anomalies with PRISM climatology
CO2	Observed ca	Constant 1901–1930	Constant 1901–1930
CO2 & T	Observed ca	CRU anomalies with PRISM climatology	Constant 1901–1930
CO2 & VPD	Observed ca	Constant 1901–1930	PRISM anomalies with PRISM climatology

Observed ca was based on the Mauna Loa annual average atmospheric CO2 record (Keeling and Whorf 2004), CRU temperature anomalies were based on the CRU-TS4.00 data product at a 0.5° spatial resolution (Harris et al. 2014). PRISM VPD anomalies were calculated from the PRISM AN81m data product at a 4 km resolution. PRISM climatology for temperature and VPD was based on the PRISM Norm81m data product at an 800 m resolution (Daly et al. 2008, 2015). Constant boundary conditions for temperature and VPD reflect the ranges and variability over period 1901–1930, without trend.

Statistical analyses

Relationships between annual climate variables and tree ring proxy were determined by correlation analyses, using the Pearson correlation coefficient (r). The independent and interactive effects of growth locality and time of growth (comparing the 100-year periods 1750–1849 and 1901–2000) on ring width and stable isotopes were tested using a two-way ANOVA. Differences between population means were tested using two-tailed Student’s t-tests at a 5% significance level. If multiple comparisons were conducted in a post hoc analysis, the significance level was adjusted using a Bonferroni correction (Cabin and Mitchell 2000). 30-year trends in growing-season temperature, growing-season VPD, preceding September-to-August precipitation and population-averaged ring widths over the period 1901–2010 were analysed by means of Pearson correlations over a 30-year running window at a 5% significance level. Periodicity in growing-season temperature, growing-season VPD, and preceding September-to-August precipitation was quantified using the FFT algorithm in PAST (Hammer et al. 2001). The association between ca and ∆ci was quantified by bootstrapping the correlation coefficient (rboot) 10,000 times over 10 ppm ca bins. For each repetition, four observations (minimum bin size) were sampled with replacement from each bin. Interpretation of the dual-isotope model was based on the direction of coordinated changes in ∆ci and δ18O. The sign of the relationships between ∆ci, δ18O and ring widths was determined with Pearson correlation on averages calculated over 10 ppm ca bins. This averaging approach was applied to prevent bias due to high sample density in the earlier part of the records. The slope of the relationship between leaf nitrogen content and maximum carboxylation capacity and the slope between ca and ∆ci were determined using standardized Major Axis (SMA) regression in R (version 3.0.2) with the SMATR package (Warton et al. 2012). Model responses to different sets of boundary conditions were analysed for their effect size, which was expressed as a percentage change relative to a common baseline response over the early-twentieth century period 1901–1930. The uncertainty of the model output as a consequence of uncertainty in the climate boundary conditions temperature and VPD was quantified with Monte-Carlo simulations consisting of 10,000 individual repetitions of each simulated year for each site. The model temperature and VPD were sampled between their minimum and maximum value for each growing-season month using a Latin square hypercube sampling strategy (McKay et al. 1979). The Latin square hypercube sampling method selects a combination of input parameters that, given a sufficiently large number of model repetitions, provides an estimate of the full range of model output associated with the sampled range of model input parameters. As the PRISM temperature and VPD anomalies are correlated (r2 = 0.25, p < 0.001), the Latin square hypercube sampling strategy weighted the sampling of these variables based on their covariance matrix.

Results

Tree-ring chronologies and twentieth century climate

Our records of annual tree ring widths revealed clear positive twentieth century growth anomalies in the treeline populations in sites 1 A and 1B, compared to the below-treeline population in site 2 (Fig. 1a). On average, the stem radii of the trees in the treeline sites 1 A and 1B increased by 93.8 (±21.9) mm and 95.2 (±35.8) mm over the period 1901–2000, respectively (two-sample t-test p > 0.05) (Fig. 1b). The average stem radius of the trees in the below-treeline site 2 increased by 33.2 (±8.8) mm over the same time period, which is significantly less than in the treeline sites 1 A and 1B (two-sample t-test p < 0.05). Stem radii of the trees in the two treeline sites showed significant growth increases between the pre-industrial period 1750–1849 and the twentieth century (two-sample t-test p < 0.05), whereas no significant change in the average stem radius was observed for the below-treeline population over this time period. An ANOVA of these century averaged growth responses showed a significant effect of site (p < 0.001) and time period (p = 0.013), but no significant interaction effect (p = 0.09). We did not observe a significant difference in DBH between populations (Fig. 1c). Considering the comparable growth responses of the treeline populations in sites 1 A and 1B, and the distinct growth responses of the below-treeline population in site 2, our following interpretation of twentieth century climate–growth responses and analyses of cellulose stable isotopes was based on annual averaged ring widths of treeline and below-treeline populations, respectively.

Figure 1.

Twentieth century growth anomalies in bristlecone pine growing in localities at the modern treeline (sites 1 A and 1B) and approximately 200 m lower (site 2). (a) Site average annual ring width (mm) and the number of trees sampled per growth locality (sample depth). (b) Pre-industrial (1750–1849) and twentieth century (1901–2000) average radial growth (mm century−1) at the three growth localities. Different letters in panel (b) indicate significant differences (p < 0.05) between groups based on two-tailed two-sample t-tests assuming equal variance. (c) Comparison of Diameter at Breast Height (DBH) of trees in the three growth localities. A one-way anova found no significant difference between the trees growing in the below-treeline site 2 (N = 6) and the treeline sites 1 A (N = 5) and 1B (N = 6) (p = 0.21). Boxes in panels (b) and (c) denote the 25th and 75th percentiles, the whiskers denote the 5th and 95th percentiles.

The analysis of climate–growth responses was focused on 30-year trends and restricted to growing-season temperature, growing-season VPD and previous September-to-August precipitation as these climate variables are most strongly associated with growth responses of bristlecone pine in the White Mountains (Salzer et al. 2009). Comparing 20th century climate trends (Fig. 2a) with trends in population-averaged growth rates (Fig. 2b) revealed that the treeline population reduced growth between the relatively cool and wet start of the twentieth century and the early-1920s, whereas the below-treeline population increased growth rates over this time period. The below-treeline population reduced growth between the 1920s and late-1950s, which may reflect a lingering response to declining precipitation until the mid-1930s and relatively warm and (high VPD) conditions persisting until the late-1930. Moreover, the period of increased rainfall that ended the 1930s drought lasted only several years. Except for this brief wetter period in the late-1930s, relatively dry and warm conditions persisted into the late 1940s. The brief period of increased rainfall in the late 1930s and early 1940s is reflected by increased stem growth in the below-treeline population, but owing to the relatively short duration, not reflected in the 30-year trends. Growth rates at treeline increased between the 1940 and 1980s in relation to a combination of warming, atmospheric drying and rising precipitation. In this time period, precipitation was especially high in the late-1960s and between the late-1970s and mid-1980s, whereas temperature and VPD were strong variable and reached maxima in 1981 and 1984. The below-treeline population showed a minor increase in growth rates in the 1960 and 1970s. Growth rates of both populations revealed statistically significant decline after the 1980s, while growing-season temperature and VPD increased, with the exception of a cold spell in the late-1990s. Preceding September-to-August precipitation showed no clear 30-year trends after the 1980s, but interannual variability was relatively high. In summary, the typical distinction between apparent temperature-limited growth at treeline and moisture-limited growth at lower elevations that was apparent between the start of the twentieth century and late 1970s dissipated during the 1980s. Declining growth rates after the 1980s appear related to the combination of warming, atmospheric drying and variable precipitation in both populations.

Figure 2.

Trends in climate variables and growth anomalies. (a) Temporal variability and 30-year trends in growing-season vapour pressure deficit and growing-season temperature anomalies (∆VPD and ∆T, respectively), and cumulative September-to-August precipitation (Pp). VPD prior to 1961 is shown with a dashed line owing to the limited direct humidity measurements in this period (Elliott 1995). (b) Temporal variability and 30-year trends in treeline and below-treeline annual radial growth anomalies (∆RW). The sign of the 30-year trends are indicated by their respective correlation coefficient (r) based on annual values in running 30-year bins. Only trends with a significance p < 0.05 are shown. Anomalies were calculated with respect to the period 1901–1930.

Cellulose stable isotopes

Signals conveyed in cellulose δ13C from our tree ring chronologies were interpreted based on the inferred changes in population-averaged δ13C-derived ci/ca. Annually pooled 13C-Suess-corrected cellulose δ13C (δ13Ccor) records are shown in Fig. S6 (available as Supplementary Data at Tree Physiology Online) for reference. We note that decadal δ13C measurements of individual trees in the population averages showed no bias due to deviating responses of individual trees (see Fig. S2 available as Supplementary Data at Tree Physiology Online). Both populations revealed positive trends in δ13C-derived ci/ca over the period 1710–2010 (correlation r = 0.72, p < 0.001 and r = 0.27, p < 0.001, for treeline and below-treeline populations, respectively) (Fig. 3a). The c/ca-ratios were notably lower in the treeline population until the early-twentieth century, with the exception of a 20-year-long positive ci/ca anomaly that occurred in the treeline population towards the end of the nineteenth century. Towards the end of the twentieth century ci/ca-ratios became comparable between populations. We found δ13C-derived ci/ca-ratios in a range between 0.35 and 0.47 for the treeline population and between 0.38 and 0.47 for the below-treeline population. These ci/ca-ratios are relatively low compared to values typically reported for leaves (Wang, Prentice, Keenan, et al. 2017), but likely reflect that δ13C values for tree ring cellulose are generally 2–3‰ higher than those of leaves (Feng 1998). Moreover, ci/ca is down-regulated in trees growing at high elevations owing to a combination of low stomatal conductance and high photosynthetic capacity (Wang, Prentice, Davis, et al. 2017). Considering century-scale responses in δ13C-derived ci/ca (comparing the period 1750–1849 to the period 1901–2000), we observed a stronger increase in ci/ca in the treeline population than in the below-treeline population (Fig. 3b). An ANOVA of these century-scale responses showed a significant effect of elevation (p < 0.001) and time period (p < 0.001), as well as a significant interaction effect (p < 0.001).

Figure 3.

Treering cellulose stable isotopes. (a) Temporal changes in δ13C-derived ci/ca measured on annual pooled treering cellulose of the treeline and below-treeline populations, respectively. (b) Century averages of δ13C-derived ci/ca reflecting the pre-industrial period (1750–1849) and twentieth century (1901–2000). (c) Temporal changes in stable oxygen isotopes (δ18O) measured on annual pooled treering cellulose. (d) Century averages of annual δ18O values. Solid lines in panels (a) and (c) reflect 5-year running averages. The correlation coefficients (r) in panels (a) and (c) indicate the sign of the trend over the period AD 1710–2010 based on annual values. Significance levels of correlations are indicated by *** for p < 0.001, and ** for p < 0.01. Boxes in panels (b) and (d) denote the 25th and 75th percentiles, the whiskers denote the 5th and 95th percentiles. Different letters indicate significant differences (p < 0.05) between groups based on two-sample t-tests assuming equal variance.

The treeline population showed a negative trend in cellulose δ18O (correlation r = −0.37, p < 0.001), whereas the below-treeline population showed a weak positive trend (correlation r = 0.16, p < 0.01) over the period 1710–2010 (Fig. 3c). An ANOVA of the δ18O records revealed a significant interaction effect between elevation and time period (p < 0.001), but no significant main effects (p > 0.05 for both factors). Moreover, the cellulose δ18O records revealed both synchronous and asynchronous multi-decadal variability (Fig. 3c). To explore the potential imprint of local climate and circulation on periodicity in cellulose δ18O, we performed a spectral analysis on both cellulose δ18O records as well as growing-season temperature, growing-season VPD, and previous September-to-August precipitation over the period 1901–2010 (see Fig. S7 available as Supplementary Data at Tree Physiology Online). Both cellulose δ18O records revealed consistent cyclicity at time scales likely corresponding to the locally relevant El Niño–Southern Oscillation (ENSO) (~2–9 year period) and the Pacific Decadal Oscillation (PDO) (~15–40 year period) (Barlow et al. 2001, Newman et al. 2003). A specifically strong periodicity of 28 years was apparent in both δ18O records. All three climate indices showed cyclicity in the ENSO band, whereas growing-season VPD and previous September-to-August precipitation also show cyclicity in the multi-decadal PDO band, but none reflect the strong 28-year signal observed in the cellulose δ18O records.

Dual isotope model

We applied the theoretical principles outlined in the conceptual dual-isotope model to interpret the combined signals conveyed in the records of δ13C and δ18O in terms of plant-physiological responses to climate variability over the period AD 1711–2010 (Scheidegger et al. 2000, Barbour et al. 2004, Grams et al. 2007). As the δ13C records are strongly influenced by plant-physiological responses to rising ca, it is somewhat problematic to discern ci/ca changes in relation to climate variability. We, therefore, detrended the records of δ13C-derived ci/ca relative to a plant-physiological baseline constant ci/ca-response to rising ca. The baseline is set at the pre-industrial (1750–1849) average ci/ca values of 0.38 and 0.41 for the treeline and below-treeline populations, respectively (Fig. 3b). The resulting time series were expressed in terms of the deviation of the time-integrated leaf intercellular CO2 concentration () from the baseline constant ci/ca-response, termed ∆ci (see Fig. S3 available as Supplementary Data at Tree Physiology Online). Our records of ∆ci revealed that increased beyond the baseline in the treeline and below-treeline populations from the pre-industrial start of the records (275 ppm ca) to the 1980s (345 ppm ca) (Fig. 4a). The positive ∆ci response observed in the treeline population was more pronounced than the response of the below-treeline population over this time period. The above-baseline responses in are halted in both populations since the 1980s, which indicates a shift to constant ci/ca-responses.

Figure 4.

Interpretation of the dual-isotope model over the period AD 1700–2010. (a) Relationship between ca and the deviation in δ13C-derived ci from the ‘active’ (constant ci/ca) response of δ13C to ca, termed ∆ci. Data points and error-bars indicate the average and one-standard deviation confidence intervals calculated over 10 ppm ca bins. The ca intervals indicated by the shaded periods (I) and (II) indicate distinct responses of ∆ci to ca (see Fig. S3c available as Supplementary Data at Tree Physiology Online); significant positive responses of ∆ci to ca were observed in period (I), whereas a shift to non-significant or negative responses was observed in period (II). Bootstrapped correlation coefficients (rboot) and their standard deviations were calculated over 10 ppm ca bins in period (I) only. (b) Temporal relationship between cellulose δ18O and ∆ci. The black boxes near the axes extremes indicate the interpretation of the adjusted dual-isotope model in terms of changes in stomatal conductance to water vapour (gsw), transpiration (E) and vapour pressure deficit (VPD), and relative changes in stomatal conductance to CO2 (Δgsc), photosynthesis (ΔA) and photosynthetic capacity (ΔAmax). (c) Temporal relationship between annual average ring width and ∆ci. (d) Temporal relationship between annual average ring width and δ18O. Markers in panels (b–d) indicate the centre of 10 ppm ca bins of the treeline and below-treeline sites. Subsequent data points are connected by solid lines, with bin-average ca indicated for every second bin. The correlation coefficients in panels (b–d) indicate the sign of the relationship calculated across the averages of 10 ppm ca bins over the period AD 1700–2010. Significance levels of correlations are indicated by *** for p < 0.001, * for p < 0.05, and NS for p ≥ 0.05.

The adjusted dual-isotope model revealed a strong negative relationship between ∆ci and cellulose δ18O in the treeline population across the length of the records (correlation r = −0.69, p < 0.05), with an exception of the period after the 1980s (with ca rising above 345 ppm) (Fig. 4b). These responses were interpreted to reflect increases in photosynthesis rates beyond the baseline constant ci/ca-response, combined with increases in transpiration until the 1980s. The relatively constant ∆ci and increase in δ18O that occurred in the treeline population after the 1980s was interpreted to reflect a transition to more modest increases in photosynthesis rates combined with a reduction in stomatal conductance and/or transpiration rates in response to rising VPD. For the below-treeline population we observed only a weak positive trend in cellulose δ18O (cf. Fig. 3c) and, therefore, found no significant association between ∆ci and cellulose δ18O (Fig. 4b). Hence, the changes in ∆ci and cellulose δ18O that occurred in the below-treeline population after the 1980s were interpreted to reflect a generally conservative transpiration strategy and, on average, photosynthetic adjustments in line with the baseline constant ci/ca response to rising ca.

We explored the relationship between population-averaged annual ring widths and cellulose stable isotopes to reveal the potential association between growth and leaf gas exchange over the period 1710–2010. For the treeline populations, we observed a strong positive correlation between ring width and ∆ci (r = 0.79, p < 0.001) (Fig. 4c) and a strong negative correlation between ring width and δ18O (r = −0.74, p < 0.001) (Fig. 4d). The continuous increases in annual growth rates and ∆ci observed in the treeline population throughout most of the nineteenth and twentieth century show a distinct change during the 1980s when growth rates started to decline and ∆ci stabilized on a constant-ci/ca trajectory. The declining growth rates since the 1980s are associated with a shift to more conservative transpiration in the treeline population. No significant correlations were observed between ring widths and the stable isotope records of the below-treeline population over this time period.

A more detailed analysis of correlations between stem growth responses, cellulose stable isotopes and climate variables over the period 1901–2010 is presented in Table S4, available as Supplementary Data at Tree Physiology Online. This analysis revealed contrasting responses of growth and gas exchange to warming-induced atmospheric drying between the below-treeline and treeline populations. No significant annual correlations were observed between ring widths and growing-season temperature, and between ring widths and previous September-to-August precipitation relationship over the 110-year period. Positive correlations were observed between previous September-to-August precipitation and ∆ci for both populations, likely reflecting an increase in stomatal conductance under increased moisture availability. The treeline population showed a significant negative relationship between previous September-to-August precipitation and δ18O, likely reflecting increased transpiration and reduced VPD during periods with increased precipitation. The below-treeline population showed no significant correlation between previous September-to-August precipitation and δ18O.

Twentieth century leaf gas exchange

Owing to the uncertainty conveyed in the ecophysiological interpretation of the dual-isotope model (Roden and Siegwolf 2012), we applied a semi-empirical leaf gas exchange model (de Boer et al. 2016) to reconstruct changes in photosynthesis, transpiration and stomatal conductance over the period 1901–2010. We applied four different combinations of model conditions to explore the sensitivity in gas exchange to changes in ca, temperature and VPD (Table 1). In the control ‘CO2, T & VPD’ simulation, modelled photosynthesis rates were comparable between treeline and below-treeline populations in the early-twentieth century (Fig. 5a). Photosynthesis rates increased consistently throughout the twentieth century in both populations, but most prominent changes were modelled for the treeline population. Our sensitivity simulations with different combinations of boundary conditions indicated that rising ca drove most of the changes in photosynthesis, as illustrated by the comparable changes observed across the four different simulations (Fig. 5b). A marginally positive effect of temperature on photosynthesis can be observed when comparing the ‘CO2’ and ‘CO2 & T’ simulations. The positive temperature effect results from average growing-season temperatures remaining below the photosynthetic optimum. The VPD boundary conditions do not affect photosynthesis in our model framework.

Figure 5.

Modelled responses of leaf gas exchange to twentieth century changes in atmospheric CO2, growing-season temperature and VPD. (a) Annual-average leaf-level time-integrated photosynthesis () as modelled in the ‘CO2, T and VPD’ simulation. (b) Modelled changes in expressed as percentage (%) with respect to the period 1901–30. The ‘CO2’, ‘CO2 & T’, ‘CO2 & VPD’ and ‘CO2, T & VPD’ simulations used different combinations of ca, temperature and VPD boundary conditions that represent twentieth century changes (see Table 1). (c) Annual-average leaf-level transpiration () as modelled in the ‘CO2, T and VPD’ simulation. (d) Modelled changes in expressed as % with respect to the period 1901–30. (e) Annual-average stomatal conductance to water vapour (gsw) as modelled in the baseline ‘CO2, T and VPD’ simulation. (f) Modelled changes in gsw expressed as % with respect to the period 1901–30. Solid lines in panels (a), (c) and (e) reflect the average of 10,000 Monte-Carlo type simulations using PRISM-based ranges in growing-season temperature and VPD as an input. Boxes in panels (b), (d) and (f) denote the 25th and 75th percentiles, the whiskers denote the 5th and 95th percentiles of the 10,000 Monte-Carlo type simulations using ranges in growing-season temperature and VPD based on PRISM as the input. Simulation medians are indicated by the circles in each box.

Modelled transpiration rates were lower in the treeline population than in the below-treeline population in the early-twentieth century (Fig. 5c). Treeline transpiration rates increased strongly between the 1910s and the mid-1930s, whereas the below-treeline population showed less pronounced variability in transpiration of this time period. Modelled transpiration rates declined between the mid-1930s and the mid-1950s and subsequently increased until the early-1970s. A transpiration minimum was modelled during the mid-1980s for both populations. Transpiration increased between the mid-1980s and 2010, with most prominent increased modelled for the treeline population. These changes in modelled transpiration rates reflect combined changes in VPD model boundary conditions (Fig. 5d) and inferred changes in stomatal conductance (Fig. 5e and f). Owing to the limited direct meteorological observations of VPD dating prior to the 1960s, VPD model boundary conditions primarily reflect the influence of temperature on VPD for the first half of the twentieth century (Elliott 1995, Seager et al. 2015). Hence, modelled changes in transpiration and stomatal conductance prior to the 1960s primarily reflect the influence of temperature on VPD.

In summary, our leaf gas exchange model showed that the altitudinally separated bristlecone pine populations responded differently to twentieth century warming and atmospheric drying. Rising ca stimulated photosynthesis in both populations, with most notable increases modelled for the treeline population. Most notable changes in transpiration and stomatal conductance were modelled for the period after the 1980s when the below-treeline population reduced stomatal conductance and the treeline population maintained a more constant stomatal conductance. Transpiration rates increased in both populations over this time period due to rising VPD and temperature. However, the decline in stomatal conductance limited the influence of rising VPD in below-treeline transpiration after the 1980s. These responses generally support our ecophysiological interpretations of the dual-isotope model and indicate that growth declines after the mid-1980s occurred in relation to a shift towards more conservative water use.

Discussion

Our absolutely dated annual chronologies of ring widths and cellulose stable isotopes reveal crucial differences the timing and magnitude of climate-induced stem growth and leaf gas exchange in bristlecone pine over the period AD 1710–2010. Our results thereby highlight that stem growth and leaf gas exchange not necessarily respond similarly to the combination of environmental changes that occurred during the twentieth century. Gas exchange specifically reflects the influence of rising CO2 on photosynthesis and rising VPD on transpiration, whereas stem growth reflects limitations associated with insufficient moisture and too low temperatures. The altitudinally separated bristlecone pine populations thereby corroborate previous work on nearby populations that reflect contrasting growth responses between treeline and below-treeline populations (LaMarche and Stockton 1974, Hughes and Funkhouser 2003, Salzer et al. 2009). The observed growth responses revealed waning temperature limitation and the onset of moisture limitation near the treeline around the mid-1980s. The shift to moisture limitation at treeline occurred approximately a decade earlier than the shift observed by Salzer et al. (2014) in a nearby site. A detailed analysis of topoclimate in relation to bristlecone pine growth by Tran et al. (2017) suggests a temperature threshold for the onset of moisture sensitivity ranging from 7.4 °C to 8 °C growing-season mean temperature. Despite considerable interannual variability in our model boundary conditions, this threshold matches the decadal-average growing-season temperature for our treeline site during the mid-1980s. Considering the comparable DBH of the populations growing at different elevations, we argue that the observed growth responses are not likely to be confounded by juvenile effects and differences in tree age. Another confounding factor could involve altitudinal differences in nitrogen inputs from air pollution (Greaver et al. 2012). However, recent evidence suggests that nitrogen loading is not significantly different between the localities sampled in the present study (Clisby 2015). Although the lack of long-term onsite meteorological measurements complicates climatological interpretation of observed growth responses, the comparable responses across bristlecone pine populations in the larger Great Basin region suggest that differential growth responses between altitudinally separated populations are most likely related to topoclimate and ongoing climate change (Tran et al. 2017).

Following the theoretical principles of the dual-isotope model (Scheidegger et al. 2000, Barbour et al. 2004, Grams et al. 2007), we inferred that the treeline population increased stomatal conductance and leaf gas exchange coeval with the positive growth anomalies between the start of the twentieth century and the 1980s. Although we interpreted the trends in stable isotope chronologies in terms of leaf gas exchange, we also observed multi-decadal periodicity in cellulose δ18O which may also reflect the influence of local climate on the isotopic signature of the soil water or changes in source water composition due to altered circulation (Gessler et al. 2014). The observed dominant period of 28 years closely matches the period of the PDO (Barlow et al. 2001, Newman et al. 2003), which is know to influence atmospheric circulation patterns over the southwester USA (Sheppard et al. 2002). Similar multi-decadal variability is not observed in the temperature, precipitation and VPD climate variables. Furthermore, our semi-empirical leaf gas exchange model supports the inferences from the dual-isotope model and highlights that transpiration and stomatal conductance increased in the treeline population during this time period. Crucially, the inferred increases in stomatal conductance occurred while growing-season VPD and ca increased. All else being equal, plants generally respond to rising VPD and rising ca by reducing stomatal conductance to prevent excess transpiration and optimize water use efficiency (Franks et al., 2013, Martins et al., 2016). Hence, the observed increases in stomatal conductance and transpiration likely reflect diminishing drought stress and improvements in tree-water status (Choat et al. 2018). The combination of warming and increased precipitation that occurred prior to the 1980s may, therefore, not only have diminished the temperature limitation to xylogenesis but also allowed the trees to increase their hydraulic transport capacity, root water uptake and leaf gas exchange, which may have further stimulated growth at the treeline.

The observed declining growth rates that occurred in the below-treeline population during periods with rising growing-season temperature and VPD suggest that, in addition to insufficient precipitation, also excess transpiration constrains growth of moisture-limited bristlecone pine. A meta analysis has shown that relatively minor increases in VPD can trigger stomatal closure to prevent hydraulic failure from excess transpiration (Bartlett et al. 2016), especially in gymnosperm species owing to their generically conservative hydraulic strategies (Meinzer et al. 2009, Carnicer et al. 2013, Manzoni et al. 2013). As suggested by Leavitt and Long (1992), the close association between temperature, VPD and evapotranspiration implies that the climate–growth responses of bristlecone pine can be confounded by warming-induced atmospheric drying. Moreover, the lagged response of the below-treeline population to increased precipitation after the 1930s drought suggests that prolonged dry periods may have a long-term imprint in bristlecone pine growth patterns, which could be related to (partial) hydraulic failure and a depletion of carbon stocks (Adams et al. 2017, Choat et al. 2018).

Our results further indicate that rising ca stimulates leaf-level photosynthesis regardless of growth locality. This generic leaf-level response is consistent with theory (Caemmerer and Farquhar 1981), experiments (Handa et al. 2005) and global modelling (Leuzinger et al. 2013). The declining growth rates after the 1980s furthermore highlight that photosynthesis is not specifically limiting stem growth in our sample of bristlecone pine (Körner 2015). Rather, our results show that twentieth century growth rates are unprecedented in the preceding three centuries, even after the 1980s. Similar exceptional growth rates were observed by LaMarche and Stockton (1974) and Salzer et al. (2009). Fatichi et al. (2016) showed that carbon fertilization can be an indirect effect of water savings resulting from reduced stomatal conductance. From this perspective, the declining stem growth rates after the 1980s suggest that any CO2-mediated water savings are relatively small compared to the influence of atmospheric drying on tree-water status. However, the changes in stem growth provide only a partial insight in tree growth responses, especially considering potential changes in root/shoot ratios (Nowak et al. 2004, Dawes et al. 2010).

We, therefore, argue that the growth and gas exchange responses observed in our sample of bristlecone pine reflect the onset of hydraulic constraints at treeline in relation to warming-induced atmospheric drying since the 1980s. Although many conifers exhibit conservative transpiration strategies under drought (Johnson et al. 2012), the extreme growth environment and specific physiology of bristlecone pine may increase sensitivity to changes in atmospheric moisture demand (Brodribb et al. 2014). Short growing seasons limit xylem development (Rossi et al. 2007) whereas ice formation in xylem sap can create embolisms that further constraint tree-water transport capacity (Pittermann and Sperry 2006). Moreover, the positive growth responses of needle lengths to summer temperatures as well as their extreme longevity likely limit the flexibility of bristlecone pine to reduce evaporative area during periods with high atmospheric moisture demand and limited water supply (LaMarche and Stockton 1974, Ewers 1982). Yet it remains uncertain whether the observed ecophysiological responses to recent climate change are indicative of relatively mild drought stress, or reflect more severe drought stress with potentially irreversible damage as a result (Choat et al. 2018). Microscopic analyses of tracheid features such as lumen and wall thickness could potentially help to interpret the severity of recent drought stress relative to responses that occurred during historic drought events (Sperry et al. 2006, Martin-Benito et al. 2013, Ziaco, Biondi, and Heinrich 2016, Rathgeber 2017).

Conclusion

The long-lived bristlecone pine is an iconic species that is highly relevant for reconstructing historic climate owing to its particular sensitivity to changes in moisture availability and temperature. Our absolutely dated annual ring-width chronologies from bristlecone pine populations growing at the treeline (~3500 m) and ~200 m below in the White Mountains of California (USA) provide site-specific information on tree growth conditions over the period 1710–2010. Measurements of stable carbon and oxygen isotopes measured on tree ring cellulose add information on changes in leaf-level gas exchange that occurred in relation to the changes in tree growth and climate. Ring widths show positive growth anomalies at treeline and consistently slower growth below treeline in relation to twentieth century warming and associated atmospheric drying until the 1980s. Growth rates of both populations declined during and after the 1980s. The cellulose stable isotopes suggest that positive treeline growth anomalies prior to the 1980s were related to increased stomatal conductance, leaf-level transpiration and photosynthesis. Reduced stem growth since the 1980s occurred with a shift to more conservative transpiration in both the treeline and below-treeline populations, whereas photosynthesis continued to increase in response to rising atmospheric CO2. Our results highlight that warming-induced atmospheric drying confounds positive growth responses of bristlecone pine populations growing near their lower temperature threshold and aggravates drought stress in populations growing near their lower moisture threshold. In addition, our results illustrate how climate change may induce relatively sudden and contrasting ecophysiological responses in an iconic conifer species.

Click here for additional data file.

Click here for additional data file.