Cluster analysis and topoclimate modeling to examine bristlecone pine tree-ring growth signals in the Great Basin, USA

Tree rings have long been used to make inferences about the environmental factors that influence tree growth. Great Basin bristlecone pine is a long-lived species and valuable dendroclimatic resource, but often with mixed growth signals; in many cases, not all trees at one location are limited by the same environmental variable. Past work has identified an elevational threshold below the upper treeline above which trees are limited by temperature, and below which trees tend to be moisture limited. This study identifies a similar threshold in terms of temperature instead of elevation through fine-scale topoclimatic modeling, which uses a suite of topographic and temperature-sensor data to predict temperatures across landscapes. We sampled trees near the upper limit of growth at four high-elevation locations in the Great Basin region, USA, and used cluster analysis to find dual-signal patterns in radial growth. We observed dual-signal patterns in ring widths at two of those sites, with the signals mimicking temperature and precipitation patterns. Trees in temperature-sensitive clusters grew in colder areas, while moisture-sensitive cluster trees grew in warmer areas. We found thresholds between temperature- and moisture-sensitivity ranging from 7.4°C to 8°C growing season mean temperature. Our findings allow for a better physiological understanding of bristlecone pine growth, and seek to improve the accuracy of climate reconstructions.


Introduction
Tree rings often reflect variability in the environmental conditions that drive tree growth. Dendroclimatology frequently uses variations in radial growth to infer interannual and longer-term variations in climate conditions (D'Arrigo et al 2001, Fritts 1966, LaMarche 1974a, Liu et al 2009. In this way, dendrochronology can be useful to help understand past climate variability, and permits inferences about climate that would be impossible with the instrumental climate record alone.
Great Basin bristlecone pine (Pinus longaeva D. K. Bailey) provides crucial long-term tree-ring data to help understand paleoclimate (Leavitt 1994, Salzer andHughes 2007), but the growth signals can be difficult to interpret. While bristlecone pine trees sometimes reach ages of nearly 5000 years (Currey 1965, LaMarche 1969, the bristlecone record presents a potentially confounding narrative because of its mixed growth signals; in many cases, not all trees at one location are limited by the same environmental variable, with some trees responding more strongly to temperature, and some being limited more acutely by moisture availability (Bunn et al 2011, Salzer et al 2009. In the semiarid high-elevation regions where bristlecone pine grows, the forest is bounded on two sides by upper and lower treelines, defining areas of the trees' dual-signal patterns. Trees at the lower forest border tend to be moisture sensitive (e.g. Hughes and Funkhouser 1998, LaMarche 1974b; tree-ring chronologies at the upper forest border have typically been used as temperature proxies (e.g. LaMarche 1974a, LaMarche and Stockton 1974, Salzer et al 2009, though topography can complicate these patterns. Salzer et al (2009Salzer et al ( , 2014 found an elevational threshold between different growth signals in bristlecone pine only 60 to 80 meters below upper treeline, with trees above that threshold primarily responding to temperature, and trees below the threshold responding to moisture availability. Further, Salzer et al (2009) conclude that only trees above that threshold show positive growth responses to recent increased warming. Trees primarily limited by temperature often differ from primarily moisturelimited trees in terms of topoclimatic conditions, the internal growth statistics of their ring widths, and the dominant frequency of variation in their ring-width growth signals (LaMarche 1974b).
Topographic variation has a significant effect on climate at a fine scale (tens to hundreds of meters), and in turn, on tree growth; we henceforth refer to climatic patterns driven by topography as topoclimate (Adams et al 2014, Bunn et al 2011, Salzer et al 2014, Slavich et al 2014. Many studies have examined topoclimate in mountain environments (e.g. Ashcroft et al 2012, Dobrowski et al 2009, Geiger 1965, Slavich et al 2014. A better understanding of topoclimate can help explain multiple growth signals in trees at a single site. Averaging ring-width measurement data from all trees in a single stand into a mean site chronology, as has been common in dendroclimatology historically (Fritts 1974, can muddle the presence of multiple growth signals in topographically complex environments.
Instead of using site chronologies to study bristlecone pine growth, the dual-signal (temperature-limited growth versus moisture-limited growth) might be explored more clearly with smaller, topographically-modified chronologies that experience different topoclimatic conditions. Additionally, because individual trees at a single site (on the order of tens of meters apart) can be limited by different factors due to topoclimate (Bunn et al 2005, 2011, Salzer et al 2014, the use of elevational lapse rates to estimate temperature might not be accurate in areas of complex terrain. At four locations in the Great Basin region, USA, we use cluster analysis to identify mixed-signa patterns in bristlecone pine radial growth, mostly using trees at the upper end of the elevational growth gradient. At sites where dual temperature-versus-moisture signals are present, we examine the different growth clusters using fine-scale topoclimate modeling and by comparing ring-width data and regional climate data. We use estimates of temperatures at tree locations to understand the fine-scale spatial differences in climate response. Our study improves upon work by Bunn et al (2011) by modeling temperatures driven by topographic variation (termed 'topoclimate'), instead of using unitless topographic indices to examine tree growth. Additionally, instead of presenting thresholds between trees of differing limiting factors in terms of elevation as Salzer et al (2014) did, we estimate thresholds in terms of topoclimate variables. By examining tree growth signals with topoclimate variables, we hope to improve understanding of physiological processes that drive and limit growth, and ultimately improve climate reconstructions using bristlecone pine.

Study areas
Tree cores and cross sections from both live trees and dead wood were collected during the summer of 2014. We attempted to choose trees across the range of complex terrain. GPS points were recorded at each tree location, and the coordinates were differentially corrected, yielding a horizontal accuracy 2 meters. The samples were added to an extant data set at four sites (figure 1) PRL,40.23 N,115.54 W, max elevation approx. 3300 m.a.s.l.). These four sites were chosen both for reasons of accessibility and to maintain continuity with past research focusing on climate-growth relationships of bristlecone pine in the same locations (Bunn et al 2011, Salzer et al 2009. Bristlecone pine grows at MWA, WHT, and PRL, and foxtail pine (Pinus balfouriana Grev. et Balf.) grows at CSL; the two species tend to have similar growth responses to climate (Bunn et al 2005).

Chronology construction
Bristlecone samples were prepared according to standard dendrochronological protocol (Stokes and Smiley 1968). This included air drying, gluing to mounting boards, and sanding with progressively finer grit. The samples were crossdated and absolute dates (with annual resolution) were assigned to the tree rings. Accuracy of the crossdating was confirmed using COFECHA software (Holmes 1999) and dplR (Bunn 2008(Bunn , 2010. The ring widths were then measured to the nearest 0.001 mm for each calendar year of the chronologies. In cases where multiple cores were collected for single trees, ring widths were averaged by tree. Though in some cases the data extend back in time to around 2000 BC, we only used data from the past four to five centuries in our cluster analyses in order to preserve robustness of both time period and chronology. For the cluster analysis, the ring-width data were detrended using a smoothing spline with Environ. Res. Lett. 12 (2017) 014007 a wavelength 2 / 3 the series length to remove non-climatic variability (Cook and Peters 1981). Because the samples in cluster analysis are exactly the same length this means each tree is detrended identically in terms of the rigidity of the spline. For the climate-growth analysis, the data extended back further in time to 1400. For those analyses, we opted to standardize the series with a conservative nonlinear model of biological growth with the form y t ¼ ae bt þ k, where a, b, and k values vary by series (Fritts 1966). We used a different detrending method for the climate-growth analysis because the data extend back further and therefore are more likely to include periods of fast juvenile growth.

Comparisons of growth clusters
Hierarchical cluster analysis was used to identify different growth signals in the tree-ring data, with growth years as variables in the cluster analysis. We used Euclidean distance matrices and Ward's method (Ward 1963) to maintain continuity with past research (see Bunn et al 2011), though other methods, such as clustering on principal components analysis, were also performed with similar results. The time periods for cluster analysis differed slightly between the four sites. The clustering method we used cannot operate with missing data; samples and time periods were chosen using the 'common.interval' function in dplR (Bunn 2008), which maximizes both the number of trees and number of years used in the analysis. We We henceforth refer to two clusters at each site as the 'red cluster' and 'blue cluster. ' We used multiple metrics to assess cluster structure (table 1), apart from visually examining dendrograms. We calculated Pearson's correlation coefficient (henceforth, 'actual r') between a mean chronology of all trees in red clusters and a mean chronology of trees in blue clusters at each site, as clusters that are more distinct from each other should be less strongly-correlated.
Using the 'sample' function in R, we randomly assigned trees to a cluster at each site, with the same number of trees per randomly-assigned cluster as in the actual clusters at each site. We then calculated a mean value chronology for each random-assigned cluster. Then, we calculated correlations (henceforth, 'random r') between those two mean value chronologies (i.e. correlation between randomly-assigned blue cluster mean chronology and randomly-assigned red cluster mean chronology) to compare to the actual  Table 1. Cluster distances and correlation coefficients. Distance and distance difference are measures of cluster quality in the analysis, with higher numbers denoting stronger two-cluster separation. A larger difference between actual r and random r suggests a stronger two-cluster separation. red-blue correlations. We performed 1000 random trials at each site, and averaged correlation coefficients from these 1000 trials to find 'random r' for each site.
Euclidean distance values at which data agglomerated into one cluster ðk ¼ 1Þ and a 'distance difference' measure, or the difference between distances at k ¼ 1 and k ¼ 2, were also calculated. These distance measures are indicative of similarity of trees in each cluster, with larger distances or distance differences denoting individual trees in a cluster being more similar to each other (i.e. better cluster 'quality').
We performed a climate-growth analysis by correlating ring-width index (RWI) for each growth cluster with regional reconstructed and modeled climate data over the period . RWI data were compared to reconstructed PDSI data from the North American Drought Atlas (Cook and Krusic 2004) and temperatures from the ECHO-G ERIK2 model (Legutke andVoss 1999, Stevens et al 2008). We gauge the climate sensitivity of these trees using climate models instead of long-term climate observations, because the models extend further back in time than the instrumental climate record, allowing us more time in common between the climate estimates and the tree-ring data. The PDSI data were reconstructed from tree-ring data, although none of the trees used in this study were used in the PDSI reconstructions. As far as we know, no chronologies from these particular sites are included in the PDSI estimates. However, it is possible there is some minimal overlap with lower elevation sites from close mountain ranges. Thus, our tree-ring data may not be completely independent of the PDSI estimates, as trees from the same region were likely used in both. Both annual and decadal correlations between RWI and regional climate data were calculated by cluster. We used wavelet analysis to examine tree growth in the frequency domain (Torrence and Compo 1998).
At sites with distinct dual-signal patterns, we compared characteristics of the two clusters in terms of topoclimate and and internal growth statistics. We modeled topoclimate using a network of temperature sensors and topographic data according to the methods described by Bruening et al (2017). This temperature sensor network included 50 sensors at each of MWA, WHT, CSL, and PRL, and collected data during 2013-2014. Bruening et al (2017) also used topographic variables to model topoclimate, including elevation, slope, eastness, southness, solar radiation loads, and topographic indices. Modeled after parameters presented by Paulsen and Korner (2014), topoclimate variables included length of growing season in days (LGS) and seasonal mean temperature in°C (SMT).
LGS is a sum of days with daily mean temperature above 0.9°C, and SMT is an average of the daily mean temperatures for the days included in the LGS (Paulsen and Korner 2014). We also used degree-hours above 5°C (DH5C) as a measure of cumulative heat sum, a common threshold for cambium cell formation in upper treeline trees (Korner 2012).
We created classification models using the topoclimate variables to predict cluster membership of individual trees. We used recursive partitioning for classification with the R package 'rpart' (Therneau et al 2014) to create these models. The models result in a single classification tree for each site, a decision tree with leaves representing the two growth signals, and branches representing climate variables that most strongly drive growth patterns (Loh 2011). At the branches, the classification trees display thresholds between the two clusters in terms of the climate variables of importance. We also compared growth clusters in terms of climate variables deemed important by the classification models with boxplots of climate variables and internal growth statistics.
All analysis was done in the R programming environment (R Core Team 2016).

Two-cluster distinction
Cluster analyses of ring-width data at MWA and WHT showed two clearly distinct growth clusters of trees (figure 2). CSL displayed only one growth signal and PRL had a more complicated structure than a dualsignal pattern (supplemental materials). Distinct clusters in the provided dendrograms are colored red and blue for easier distinction and for reference as 'red cluster' and 'blue cluster.' Red and blue clusters at MWA and WHT were lessstrongly correlated with each other than clusters at CSL and PRL (table 1). The correlation coefficients (random r) between clusters when trees were randomly assigned were significantly higher than the actual red-blue correlation coefficients at MWA and WHT (table 1). On the other hand, random r values were similar to actual r values at CSL and PRL. Note that our reported correlation coefficients are not correlations between individual trees within a cluster, but are rather correlations between mean chronologies of different clusters. Cluster analysis at MWA and WHT each produced more distinct clusters in terms of distance and distance difference than CSL and PRL. Therefore, we continued with analysis only at MWA and WHT, the two sites that showed clear dual-signal patterns.

Climate-growth analysis by cluster
Correlations between ring-width index (RWI) and climate data both as annual data and with 20-year smoothing splines of those data are displayed in table 2. At both sites, relationships between 20-year splines of RWI and temperature have the highest coefficients for the blue clusters. For red clusters, there are notable correlations between RWI and annual PDSI data.
Environ. Res. Lett. 12 (2017) 014007 Wavelet analysis plots (figures 3 and 4) show that variation is most apparent at multidecadal frequencies in the blue cluster signals, and more sub-decadal variation is present in red cluster signals. Additionally, blue clusters at both sites show increased growth in recent years, as displayed by increases in RWI.

Cluster topoclimate and growth differences
Trees at MWA and WHT show distinct two-cluster patterns in radial growth. We used several methods to compare cluster characteristics at the two sites: boxplots of topoclimate variables and growth characteristics were useful in comparing clusters, and we created classification models to find the most important topoclimate variables in determining cluster membership.
At both MWA and WHT, trees in the blue clusters grow in areas with colder seasonal mean temperature (SMT), shorter growing seasons, and fewer annual degree hours above 5°C (figure 5). Blue cluster trees tend to have higher first order autocorrelation and higher standard deviation of RWI.
The second method of comparing growth signals/ clusters by topoclimate used R package 'rpart' to create models that predict tree cluster membership based on topoclimate variables. Classification trees in figure 6 display the most influential topoclimate variable in determining cluster membership at both sites. At MWA, trees growing in locations with SMT < 8 C tend to be in the blue cluster, while trees in areas with SMT > 8 C tend to be in the red cluster. At WHT, a similar pattern is evident, but with a threshold of SMT ¼ 7:4 C instead of 8 C. Cohen's kappa statistic was used as a measure of the quality of the classifications (Cohen 1960). The kappa value at MWA ðk ¼ 0:67Þ  Figure 2. Dendrograms of clustering from MWA (above) and WHT (below). There are two distinct growth signals at each site, henceforth the red and blue clusters. Table 2. Correlations between RWI of clusters and climate variables from climate-growth analysis over 1400-1990. PDSI estimates are reconstructed from tree ring data, and temperature estimates come from the ECHO-G ERIK2 model. Most notable coefficients are in bold. MWA is above, WHT is below. We did not perform significance testing, because the large sample sizes would likely result in finding meaningless statistical significance (Sullivan and Feinn 2012   Most variability is low-frequency (decadal to multidecadal) in the blue cluster and high-frequency (annual to sub-decadal) in the red cluster.
Environ. Res. Lett. 12 (2017) 014007 was quite high, suggesting substantial agreement between predicted class and actual class in the classification models. The kappa value at WHT ðk ¼ 0:52Þ was lower, suggesting moderate agreement between predicted and actual classes. Univariate classification models were also created, using only one topoclimate variable as a predictor at a time. The results of these univariate models are displayed in table 3; thresholds between clusters for the three topoclimate variables are quite similar at the two sites. The spatial distribution of red cluster and blue cluster trees is displayed in figure 7. The SMT thresholds of 8 C and 7:4 C at MWA and WHT, respectively, are shown as black contour lines in figure 7, and distribution of trees dictated by these thresholds is evident.

Two-cluster distinction
The presence of multiple growth signals at a site is potentially significant, as it would justify the need for analysis based on growth signal (likely related to limiting factor) rather than overall site chronologies in dendroclimatological studies (Fritts 1976). When reconstructing past temperatures with tree rings, the assumption that all trees at a site are temperaturelimited causes the introduction of noise from the moisture-limited trees that might get mixed in, thus diluting and masking the temperature signal in the chronology.
The lack of a clear two-signal growth structure at CSL and PRL could be a result of sampling techniques when collecting wood or a true one-signal growth pattern at the site. Samples were collected over a much narrower range of temperatures (and elevations) at  CSL and PRL than at MWA and WHT. The lack of diversity of sampling locations likely contributed to a lack of a dual-signal structure at CSL and PRL.

Climate-growth analysis by cluster
Several lines of evidence suggest that low-frequency signals are associated with temperature variability and that high-frequency signals are related to moisture variability in bristlecone chronologies from the same site (LaMarche 1974b, Esper et al 2002, Hughes and Funkhouser 2003. Low-frequency temperature signals are often apparent on multidecadal and multicentennial terms, while high-frequency moisture signals tend to be observed with one to several year variations (Hughes and Funkhouser 2003). Blue clusters at MWA and WHT show most power around 32-years and greater periods. Red clusters show higher power in higher frequencies, though there are some low-frequency signals mixed in (figures 3 and 4). Frequency analysis (figures 3 and 4) and correlations between RWI and regional climate data (table 2) suggest that blue clusters at both sites are made up of temperature-sensitive trees, and red clusters at both sites are composed of moisture-sensitive trees. We will henceforth refer to the blue and red clusters as temperature-sensitive and moisture-sensitive clusters, respectively.
At both sites, the temperature-sensitive clusters are accurately composed almost exclusively of temperature-limited trees. Temperature data from the ECHO-G ERIK2 model correlate very well with RWI with 20-year smoothing splines (r = 0.66 and r = 0.74 at MWA and WHT, respectively) for the temperaturesensitive clusters (table 2). These clusters also have negligible correlation with PDSI, suggesting that trees in the temperature-sensitive clusters were properly assigned. Because the temperature-sensitive cluster has less of a mixed signal than the overall site chronologies, mean chronologies from the temperature-sensitive clusters in this study would be appropriate for temperature reconstructions. Cluster analysis did not produce as clean of results for moisture-sensitive clusters, as there appear to be some temperature-sensitive trees within the moisturesensitive clusters at both sites. Bunn et al (2011) addressed the more complicated moisture-sensitive signal by using more than two clusters in their analysis.
We followed a conceptual model of only two limiting factors on tree growth (temperature and moisture availability) for simplicity. One of the strengths of this study is the use of fine-scale temperature data, and a definitive temperature-sensitive cluster allows use of those data; the same data are not currently available for moisture. Tree growth response to climate can change over the lifespan of a tree, both in magnitude of correlation coefficient with climate variables, and in changing from positive to negative relationship or vice versa (Jacoby and D'Arrigo 1995, Biondi 2000, Sullivan et al 2015, Zang and Biondi 2015. Nonstationary climatetree growth relationships could be especially applicable with bristlecone pine because of its extreme longevity. The use of wavelet analysis allows frequency analysis over a multicentennial time period. To avoid mixed signals, consideration of stable relationships between growth and climate should be considered prior to the application of these types of data in climate reconstructions. Paulsen and Korner (2014) use length of the growing season (LGS) and growing season mean temperature (SMT) to predict treeline position in a global model. We apply the same variables, along with annual degree hours above 5°C (DH5C) as a cumulative measure of temperatures experienced by trees. The global TREELIM model by Paulsen and Korner (2014) found that a minimum LGS of 94 days and a minimum SMTof 6.4°C were required for tree growth. A more regional model focusing on bristlecone pine by Bruening et al (2017) found that at sites in the Great Basin region, minimum LGS ranging from 147-153 days and minimum SMTof 5.5°C-7.2°C dictated upper treeline position. Instead of thresholds for growth possibility, the values presented here (figure 6 and table 3) are estimated thresholds between temperature sensitivity and moisture sensitivity, similar to those alluded to in terms of elevation by Salzer et al (2014). Naturally, these values are greater than the minimum limits for growth presented by Paulsen and Korner (2014) and Bruening et al (2017).

Cluster topoclimate and growth differences
At both MWA and WHT, the creation of classification trees found that SMT is the variable that most distinctly defines differences between temperature-and moisture-sensitive trees. It is possible to predict whether trees will be part of a temperature-sensitive cluster or a moisture-sensitive cluster based solely on SMT as a predictor variable for the majority of trees at MWA and WHT. Bunn et al (2011) characterized clusters of trees in terms of unitless topographic indices. The use of modeled temperature values in our study allows for a better physiological understanding of bristlecone pine growth response to climate than thresholds based on elevations or topographic indices allow, as temperature is the primary driver of growth in many high-elevation trees (LaMarche 1974a, Bunn et al 2011, Salzer et al 2009. Warming climates will likely further complicate the puzzle of understanding mixed signals in bristlecone pine growth. Lloyd and Graumlich (1997) suggest that with warming temperatures, precipitation will play a larger role in dictating upper treeline position and structure. Complicated topography in mountainous areas often provides climate refugia (Scherrer and Korner 2011), but topoclimatic conditions in these areas will likely impact growth patterns in bristlecone pine. Temperature thresholds between limiting factors of growth like the ones we present in our study will be important in interpreting mixed growth signals under changing conditions.
Internal tree-ring statistics such as standard deviation and first order autocorrelation coefficient (AR1) are sometimes indicative of different limiting factors of a tree-ring chronology. Ring-width data from MWA and WHT demonstrated higher AR1 for temperature-sensitive trees, and we also observed higher values of standard deviation in temperaturesensitive trees, similar to what Wilson and Luckman (2003) found with upper treeline Engelmann spruce. The separation between clusters for standard deviation and AR1 was not as defined as cluster separation for temperature variables (figure 5).

Conclusions
This research expands on work by Bunn et al (2011) and Salzer et al (2009Salzer et al ( , 2014 by using modeled topoclimatic variables instead of variables that are physiologically less-influential, such as elevation or topography to examine multiple growth signals in bristlecone pines. Cluster analysis is a strong tool for identifying patterns in multivariate data, and in this study, it was used to find groups of trees with similar limiting factors of growth. A central goal of this study was to more easily identify trees that are sensitive to a certain limiting factor to improve climate reconstructions of that factor. At two of our four sites, we found two different growth signals (clusters) within the many trees sampled. One subset of trees showed a temperature signal while the other subset of trees displayed a mode of growth consistent with a precipitation signal. Our findings largely coincided with results of past research (LaMarche and Stockton 1974, Bunn et al 2011, Salzer et al 2014: temperature-sensitive trees have the tendency to grow at higher elevations and in colder temperatures than moisture-sensitive trees. However, using topoclimate modeling, we were able to take these conclusions one step further. We were able to define thresholds between trees limited by temperature and trees limited by moisture availability by differences in growing season mean temperatures. Attempts made to use potential evapotranspiration and topographicallyderived drought indices as water balance proxies to improve understanding of hydroclimatic thresholds were unfruitful. Future studies could use hydrologic modeling to better understand moisture limitations on bristlecone pine growth. A possible candidate for an appropriate hydrologic model is the California Basin Characterization Model (Flint et al 2013).
One aim of this study was to help improve the accuracy of climate reconstructions that use bristlecone pine. Cluster analysis and topoclimate modeling can be valuable and powerful tools when used to better understand the complicated bristlecone pine growth signal. Detailed descriptions of different growth signals based on topoclimate will help researchers better target trees for sampling and more clearly understand tree physiology.