Wood samples were collected from bristlecone pines growing at 3355 m elevation in the Snake Range of eastern Nevada, within the Great Basin of North America (Figure 1A). The study site (38°54′22″N, 114°18′32″W) is part of the Nevada Climate-ecohydrological Assessment Network (NevCAN), a mountain observatory which consists of two valley-to-mountain transects established between 2010 and 2013 to collect long-term data on climate variability and its ecohydrological impacts (Mensing et al., 2013). Bristlecone pine is the dominant overstory species at the study area, with scattered individuals of limber pine (Pinus flexilis E.James) and Engelmann spruce (Picea engelmannii Parry ex Engelm.) also present. Given the short record (2011–2015) of sub-hourly climatic data recorded by the NevCAN station, monthly summaries of total precipitation, average, maximum, and minimum temperature for the site were obtained from the 4-km gridded PRISM (Parameter-elevation Relationships on Independent Slopes Model) dataset (Daly et al., 2008), and interpolated at the coordinates of the study site. During the 1895–2013 period, average annual air temperature was 3.9°C (Figure 1B) and mean total annual precipitation was 628 mm, showing an annual cycle dominated by cool-season precipitation, with summer being both the warmest and the driest period.

Ten bristlecone pines (P. longaeva D.K. Bailey) of similar size and vigor were selected. The average diameter at breast height (1.3 m from the ground) was 38.5 ± 8.6 cm while tree height was 10.2 ± 2.0 m. Two increment cores were collected at breast height from each tree using a Swedish increment borer with inner diameter of about 4.3 mm (Grissino-Mayer, 2003). In the laboratory, these 20 cores were mounted, sanded, and crossdated using standard dendrochronological procedures (Stokes and Smiley, 1996). Ring-widths were measured to the nearest 0.001 mm with a Velmex measuring system interfaced with the program MeasureJ2X. Crossdating of ring-width measurements was checked with the dplR library (Bunn, 2010) for the R numerical software environment (R Core Team, 2015).

A site ring-width chronology was developed using all 20 cores extracted from the ten selected trees (Table 1). Measurements of anatomical parameters were performed on a subset of ten increment cores (one per each tree) using a non-destructive procedure with a Confocal Laser Scanning Microscope system that captures microscopic images directly from wood surfaces (Liang et al., 2013a). While the application of CLSM methods is well established in wood sciences (Knebel and Schnepf, 1991; Singh and Donaldson, 1999), they have been used only recently to develop time series of anatomical parameters for dendroclimatic purposes. Sample preparation followed standard procedures (Liang et al., 2013a), starting with the application of a non-Newtonian fluid (i.e., corn starch) on each mounted core to fill cell lumens and prevent cell walls from collapsing during the cutting operation (Schneider and Gärtner, 2013). Each core was then surfaced using a WSL core-microtome (Gärtner and Nievergelt, 2010), and a safranin solution was applied with a brush to the core surface to increase the contrast between cell walls and lumens. Light from a helium neon laser with wavelength of 543 nm generated by an Olympus FluoView FV300 CLSM system activated auto-fluorescence of wood. The reflected light passed through a confocal aperture where it was transformed in an electric signal, and composed into a distortion-free, well contrasted digital image at 100× magnification (Liang et al., 2013a).

Successive images from a core were merged with Adobe Photoshop^TM and a shading correction was applied to remove background auto-fluorescence noise (Moëll and Donaldson, 2007). Digital CLSM produced images with black background and green cell structures (Figure 2A). Composite images were analyzed with the WinCELL^TM software (Guay, 2013). Tree-ring boundaries were visually identified on the microscopic sections and anatomical features were measured for each ring (Figure 2B). Since latewood was in general extremely narrow (i.e., one or two rows of cells), cell measurements were aggregated into a single annual value. Measurements were checked so that errors related to sample imperfections (e.g., resin ducts or surface breaks) could be excluded using either an automated procedure (i.e., by setting proper filters) or by the operator (Liang et al., 2013a). Annually resolved time series of anatomical parameters included average LA, the average LA of the 30% largest cells (LA30), lumen radial diameter (LD), lumen tangential width (LW), double cell wall thickness (DWT), and cell radial diameter (CD, which was given by the sum of LD and DWT) (Table 1). To perform dendroclimatic comparisons and considering the fact that typically many more ring-width samples can be processed in the same time required to analyze a certain number of anatomical parameters, an additional ring-width chronology was produced using only the 10 cores that were analyzed for anatomical parameters (TR10).

Annual tree-ring parameters, either width or anatomical features, were combined into a site chronology after removing age-related trends and other non-climatic variability (Cook and Kairiukstis, 1990). A cubic smoothing spline (Cook and Peters, 1981) was fit to each series of ring-width and cellular measurements, and indices were obtained as ratios between the measurements and the corresponding spline values. The median of all indices available for a year was used to produce the chronology value for that year, as follows:

where Ī_t = chronology value in year t = median annual index; n_t = number of samples in year t, with n_t ≥ 5; w = crossdated ring parameter of sample i in year t; y = value of sample i in year t computed by fitting a cubic smoothing spline with 50% frequency response at a period of 50 years; w_t/y_t = dimensionless index value of sample i in year t. Empirical measures of dendroclimatic signals (Hughes et al., 2011), in particular average inter-series correlation (RBAR_EFF), expressed population signal (EPS), and signal-to-noise ratio (SNR), were computed to test the strength of the environmental information embedded in the standardized chronologies using the maximum overlap of pairwise correlations (Bunn et al., 2014). The all-lag sensitivity of the chronologies was quantified by the Gini coefficient (Biondi and Qeadan, 2008).

Response functions between tree-ring chronologies and monthly values of average, maximum, and minimum air temperature and total precipitation from the PRISM dataset were computed using the R package treeclim (Zang and Biondi, 2015). The period 1950–2013 was selected to avoid climate data uncertainty, mostly caused by very sparse coverage during earlier years (1895–1949) in the western USA (Guttman and Quayle, 1996). Dendroclimatic relationships over the water-year window, from the previous October to the current September, were tested using 1000 bootstrapped samples. Moving response functions (Biondi, 1997) for a 30-year sliding interval allowed for evaluating the stability of dendroclimatic relationships.

Seasonal correlations between climate and tree-ring chronologies were computed by means of partial correlation for season lengths from three to six months, using the approach of (Meko et al., 2011) as implemented in treeclim (Zang and Biondi, 2015). Skill of ring-width and anatomical proxies for reconstructing past climate was tested by performing ordinary least squares regression for the most significant seasonal correlation that emerged from partial correlation analysis between ring-width and anatomical parameters against precipitation. The period 1950–2013 was used for a split calibration-verification exercise using the reduction of error (RE), the coefficient of efficiency (CE), and the Durbin–Watson test (Biondi et al., 2008).

Stem ages of sampled bristlecone pines ranged from 121 to 240 years. Crossdating was simplified by the absence of false or missing rings, and by the presence of several “pointer years” (e.g., 1966, 1924, 1902, 1884). In all cores the 1902 earlywood showed the typical features of “frost rings” (Schweingruber, 2007), indicative of a late freezing event in the spring of that year. In the 1884 latewood an early-frost event was visible in about half of the cores. Without detrending of individual series, median chronologies of ring-widths (TR20 and TR10) showed an overall declining pattern, whereas median chronologies of anatomical parameters (LA, LA30, LD, CD, LW, DWT) presented an opposite trend, in particular DWT and LW (Table 1). Most anatomical parameters were not correlated with ring-widths (either TR20 and TR10), with the only exception of LW and DWT (Table 1). Standardized chronologies started in 1870, the first year with at least five overlapping anatomical samples, and ended in 2013, for a total length of 144 years (Figures 3A–H). Anatomical chronologies presented similar or lower first-order autocorrelation (ar1) compared to ring-width series. The parameters LW and DWT again behaved differently, presenting the highest ar1 values, while CD had the lowest one (Table 2).

Common environmental signals recorded by the tree-ring chronologies decreased from ring-width to wood anatomy parameters. Measures of empirical signal strength (Table 2) for the site ring-width chronology (TR20) were in the same order of magnitude of other P. longaeva chronologies (LaMarche and Stockton, 1974). Signal strength in TR10 was also high and values of RBAR_EFF, SNR, and EPS exceeded the minimum thresholds suggested for sufficiently replicated chronologies (Wigley et al., 1984). Among cellular features, LD showed the highest RBAR_EFF, SNR, and EPS (Table 2). Parameters linked to plant’s hydraulic performance (LA and LA30) and CD (a combined parameter including LD and DWT) had lower SNR and EPS, but higher Gini coefficients than other anatomical chronologies (Table 2). Noise was predominant (SNR < 1) in chronologies of lumen width (LW) and DWT, and both RBAR_EFF and EPS dropped to extremely low values (Table 2). These two anatomical parameters were therefore excluded from further dendroclimatic calibration.

Dendroclimatic response functions highlighted the presence of a general positive correlation with current year precipitation in 1950–2013 (Figure 4). All the anatomical parameters showed a significant correlation with March and July precipitation (Figures 4C–F). Radial size of cell lumens (LD) and average LA were also positively related to August precipitation (Figures 4C,E). Only the ring-width chronology based on 20 cores (TR20) showed a significant inverse relationship with September minimum temperature (Figure 4A), while LA presented a positive correlation with previous October minimum temperature (Figure 4E).

Since no significant relationships emerged between anatomical parameters and mean, maximum, or minimum temperature (with the only exception of the above mentioned relationships between LA and previous October temperature), moving response functions were computed only for precipitation. The importance of March and July precipitation for shaping wood cellular features was particularly evident for lumen and cell diameter (CD; Figures 4C–D). March was stronger and more temporally stable in LD and CD chronologies compared to other anatomical parameters (e.g., LA). Response to July moisture remained positive from 1950 to 2013, even though it seemed to increase in LD, CD, and LA in the last 35–40 years (Figures 4C–E). Approximately in the same period, the precipitation signal in April and June reversed for both LD and CD, showing a positive effect in April and a negative effect in June.

Partial correlation analysis for the period 1950–2013 indicated a positive seasonal precipitation signal for anatomical chronologies (Figure 5). Ring-width chronologies presented a significant seasonal relationship with total precipitation over the previous autumn and winter (October–March), which emerged in both TR20 and TR10 (Figures 5A–B). Anatomical parameters LD and CD showed a higher correlation with spring-summer precipitation during the current year. Lumen diameter (LD) was highly correlated with March–August precipitation (r = 0.51; Figure 5C), while CD was positively linked to total precipitation from May to September (r = 0.46; Figure 5D). LA (Figure 5E) and the average size of the top 30% largest cells (Figure 5F) were correlated to winter and spring precipitation, LA from January to May (r = 0.46) and LA30 from previous October to current March (r = 0.45).

Lumen diameter emerged as the best potential proxy for reconstructing past seasonal precipitation from 1870 to 2013 (Figure 6; Table 3). The regression model between LD and total March–August precipitation explained about one-third of the total variance (R² = 0.31), and verification indices (both RE and CE > 0) were satisfactory for the two split calibration/verification periods (1950–1981 and 1982–2013; Table 3). Verification indices between CD and May–September total precipitation were also acceptable (in particular for the period 1982–2013) but the model built using CD explained less variance (R² = 0.23) and the Durbin–Watson test pointed to higher residual autocorrelation (Table 3). Other anatomical parameters (LA and LA30) had good verification statistics, respectively, for the period 1950–1981 and 1982–2013, but the overall model explained a limited amount of variance (R² = 0.20 and 0.19 for LA and LA30). Ring-width chronologies could not be used to reconstruct past variability of total October–March precipitation, as the verification statistics and model’s R² were extremely low for both TR20 and TR10 (Table 3).

These 144-year long, crossdated, cellular chronologies represent a first attempt at developing proxy climate records from anatomical features of Great Basin bristlecone pine. Similarly long time series of anatomical parameters (e.g., cell wall thickness; vessel area) have already been generated in Canada and Europe for hardwood (Tardif and Conciatori, 2006; Fonti and García-González, 2008) and conifer species, including Pinus nigra Arn. (Martin-Benito et al., 2013), Pinus sylvestris L. (Eilmann et al., 2009; Seo et al., 2012; Pritzkow et al., 2014), and Picea abies Karst. (Castagneri et al., 2015), Picea mariana Mill. (Wang et al., 2002), and Picea glehnii Mast. (Yasue et al., 2000). To our knowledge no such chronologies exist in the western US, despite the presence of several long-lived conifer species. In particular, the importance of P. longaeva for paleoclimatic reconstructions can hardly be overstated, given the species’ property of attaining stem ages in excess of 5000 years (Currey, 1965).

Xylem cellular features were not correlated with ring-width, with the only exception of the cell wall thickness and tangential LW. Our results then suggest that wood anatomy of bristlecone pine carries a specific environmental signal, which is not a replicate of the signal encoded in ring-width series, but rather represents a source of additional paleoclimatic information (Eckstein, 2004). In our dataset, ar1 of most anatomical parameters was lower than for ring-width chronologies, as intra-annual tracheids are largely affected by environmental conditions at the time of their formation, in particular the tree hydration status (Steppe et al., 2015). Reduced ar1 in anatomical time series (i.e., tracheid and vessel size) compared to ring-width series has been in fact reported in conifers (Pritzkow et al., 2014; Pacheco et al., 2016) and in hardwoods (García-González and Eckstein, 2003; Fonti and García-González, 2008). The values of ar1 for the chronologies of DWT and LW can be explained considering the seasonal dynamics of wood formation. Cell wall thickness can present high temporal autocorrelation (Martin-Benito et al., 2013), likely because of how previously stored carbon pools contribute to the current year growth (Hansen and Beck, 1990), while the tangential diameter of developing xylem cells remains almost constant during the elongation phase (Cuny et al., 2014).

The EPS of LD and CD was comparable with values found in ring-widths of bristlecone pine growing at the treeline in the southwest US (Salzer et al., 2014), yet lower than the common signal in ring-widths (TR10). Weaker empirical signal strength, in particular inter-series correlation, for cellular chronologies compared to ring-width series, can be the consequence of lower inter-annual variability in microscopic features of wood (Olano et al., 2012; Liang et al., 2013b; Pritzkow et al., 2014). On the other hand, a relatively weak common signal in annual anatomical parameters can correspond to significant relationships with intra-annual environmental variables (Yasue et al., 2000). Information recorded in xylem anatomy may in fact reflect internal physiological processes (Fonti and García-González, 2008) rather than limiting factors exerted over the entire growing season (Eckstein, 2004).

Climate-growth relationships at the cellular level have been investigated using undetrended series (Liang et al., 2013b) as well as standardized chronologies (Bryukhanova and Fonti, 2013). Standardization is a critical issue in dendrochronology (Hughes et al., 2011), and its importance has become evident even for xylem anatomy, especially when the length of individual series exceeds multiple decades (Pritzkow et al., 2014, 2016; Carrer et al., 2015). Our cellular chronologies of bristlecone pine were characterized by an age-related trend, which made standardization useful. Cellular parameters were better dendroclimatic proxies than ring-width chronologies when detrending techniques commonly applied to ring-width series were used. In our study, dendroclimatic response functions of annual tree-ring parameters showed a dominant precipitation signal, especially for anatomical parameters LD, CD, and LA, and a very weak temperature signal. This is in contrast with the positive response to temperature and the weak correlation with precipitation found by (Salzer et al., 2009) at the Snake Range. Such discrepancy could be related to the age of the individuals sampled in these studies, and it reflects the complexity of climatic responses near treeline (Jolly et al., 2005), which can vary even at very short spatial scales, especially for P. longaeva (Salzer et al., 2014).

The ecological meaning of climate-anatomy relationships can be better understood if integrated with information on xylem phenology (e.g., the timing of wood formation) and intra-annual dynamics of stem growth obtained by repeated microcoring and dendrometer records (Rossi et al., 2012; Cocozza et al., 2016; Pacheco et al., 2016). In a study recently conducted at the same site (Ziaco et al., 2016), production of new tracheids in bristlecone pine during two consecutive years (2013–2014) started at the beginning of June, immediately after snowmelt, peaked in late June, and wood formation ended in mid-September. Cellular measurements showed that mean air temperatures between 6°C and 7°C were linked with the production of new tracheids, while dendrometer records highlighted the influence of soil temperature for snowpack thawing and subsequent plant-soil water exchanges (Ziaco et al., 2016).

Conifer species growing in drought-prone environments escape hydraulic failure (e.g., cavitation) because their reduced conduit size (e.g., LA, LD) and enhanced cell wall thickness (DWT) can resist higher negative water pressures (Hacke et al., 2001). The positive and temporally stable precipitation signal found in cellular chronologies of LD, CD, and LA can be divided in a winter (i.e., March) and summer component (i.e., July–August), which most likely relate to water-use strategies throughout the growing season. March precipitation at the study area is important to recharge soil water, providing trees with moisture needed in the first phases of xylogenesis for maintaining cellular turgor (Turcotte et al., 2009). Summer precipitation relieves water stress and increases the margins for hydraulic safety, leading to the formation of larger cellular elements (Martin-Benito et al., 2013; Ziaco et al., 2014; Pacheco et al., 2016). It should also be noted that at high elevation (i.e., near treeline), reduced lumen sizes can also be interpreted as an adaptation of bristlecone pine to reduce the risk of freezing-induced embolisms (Pittermann and Sperry, 2003). These ecophysiological processes make it possible to use the inter-annual variability of tracheid size as a reliable descriptor for tree hydraulic adjustments (Bryukhanova and Fonti, 2013).

The seasonal climate-growth relationships uncovered by partial correlation analysis suggest that in these arid and semi-arid environments water availability plays a crucial role on the xylem anatomy of bristlecone pine. An adequate supply of moisture is needed during all phases of cellular differentiation, especially in dry environments (Vieira et al., 2014), and the anatomical structures formed throughout the growing season are often related to intra-annual climatic episodes (Deslauriers and Morin, 2005), including the most extreme ones (Carrer et al., 2016). Water availability also controls carbon consumption and modulates internal distribution of non-structural carbohydrates (Sevanto et al., 2014), ultimately driving mortality patterns of entire forest ecosystems (McDowell et al., 2008). Different time scales of physiological processes driven by water availability might therefore explain the prolonged seasonal response of xylem anatomy identified in this study on P. longaeva. Even the negative correlation between ring-width chronologies and September temperature could reflect an indirect response to moisture, as suggested by the positive correlation between TR20 and September precipitation that emerged after 1970 (Figure 4A). At high elevations, daily temperature variability is regulated by incoming solar radiation, therefore cloudy skies are associated with a higher chance of precipitation, limited evapotranspiration, and lower air temperatures (Biondi and Rossi, 2015). Water is a key driver of cambial resumption in the dry season (Camarero et al., 2010), and trees could therefore use the additional moisture available in late August-early September to delay the cessation of wood production.

Bristlecone pine, thanks to its remarkable longevity, has provided reliable reconstructions of temperature and precipitation in the western US over the last millennium, through chronologies based on ring-widths (Routson et al., 2011; Tolwinski-Ward et al., 2015) or stable isotopic ratios (Leavitt, 1994; Bale et al., 2011). In this study, ring-width and anatomical chronologies performed differently as dendroclimatic proxies for reconstructing past seasonal precipitation. Chronologies of LD, CD, and LA had empirical signal strength statistics lower than both TR20 and TR10, but provided better precipitation proxies. Our reconstruction of total March–August precipitation based on LD of tracheids passed typical calibration/verification tests, even if developed using a limited amount of samples (10 cores). This result suggests that dendroclimatic information recorded in cellular anatomy of bristlecone pine could be used to refine our understanding of past climatic variability in the Great Basin, including the frequency of intra-annual extreme events or the occurrence and duration of severe drought episodes.