From treeline to species line: Thermal patterns and growth relationships across the krummholz zone of whitebark pine, Sierra Nevada, California, USA

ABSTRACT

A multiyear study of forest-to-alpine ecotones across extensive krummholz zones in whitebark pine (WBP; Pinus albicaulis), Sierra Nevada, California, resolved mean treeline growing season temperature (GST) of 9.3°C, 2.6°C warmer than global thresholds previously described, and mean growing season length of 143 days. Temperatures declined with increasing elevation; GST at the upper krummholz line (8.9°C), however, was 2.2°C warmer than the mean global treeline threshold, suggesting that by thermal criteria these environments should support tree growth. Possible explanations for the warm conditions and persistence of krummholz rather than treeline advance include a role for moisture limitations, disequilibrium with Little Ice Age temperatures, and the influence of krummholz as a buffer to treeline dynamics. Radial growth in treeline WBP trees was negatively correlated with maximum annual temperature and positively correlated to water year precipitation. Krummholz stems had low correlations to one another, to treeline trees, and to climate, suggesting nonclimatic controls on growth. These findings underscore the variable nature of treeline response to climate change, suggest that krummholz ecotones behave differently from diffuse treelines, and add to examples of mountain conifers that may be exhibiting lag effects and have not shifted with contemporary warming.

KEYWORDS:

• Whitebark pine
• Pinus albicaulis
• treeline
• krummholz
• Sierra Nevada
• radial growth

Introduction

Mountain ecosystems are widely assumed to have characteristic vulnerabilities to changing climates (Beniston 2003; IPCC 2019). Isotherms are expected to shift up as temperatures rise, with concomitant upward migration of species as they seek favorable climates (Bell, Bradford, and Lauenroth 2014). For tree species, the alpine ecotone has long been recognized as responsive to climate, although controls are complex and continue to elude comprehensive understanding (Malanson et al. 2011, 2019; Camarero et al. 2017). The nature of treeline response has been informed by global studies that elucidated factors consistently correlating with limits of upper elevation tree growth (Körner and Paulsen 2004; Paulsen and Körner 2014). Körner and colleagues emphasized the role of upright tree stature at treeline for its coupling of the crown with free-air temperatures. Climatic treeline, where atmospheric conditions enforce the limits of tree occurrence, occurs at the highest extent of groups of trees growing ≥3 m in height and growing season ≥94 days (Körner and Paulsen 2004). Under those conditions global treeline coincided with mean growing season temperature (GST) of 6.7°C. Above climatic treeline meteorological conditions change abruptly as a result of upright tree absence. Low-stature vegetation becomes decoupled from regional air temperatures and microclimates increase in importance (Grace 1988; Holtmeier 2009). Plants in the alpine zone survive harsh conditions through multiple adaptations, with many tree species becoming stunted and shrubby, a form known as krummholz (“crooked wood”). The upper limits of species occurrence defines the upper species line (Körner 2012a). For conifers that form krummholz, species line can occur far above climatic treeline.

Although GST has been demonstrated to control many treeline limits, other factors influence treeline dynamics. A global review of 166 treeline sites suggested that diffuse treeline ecotones, where upright trees diminish in size as they approach their high-elevation limit, were most sensitive to GST and responsive to contemporary warming (Harsch et al. 2009). For example, undisturbed and diffuse treelines of Scots pine (Pinus sylvestris) showed rapid increase in density and movement upslope during the latter half of the twentieth century in response to summer warming (Kullman 2005). Krummholz-type treeline ecotones, by contrast, showed greater stability in the face of warming and appeared to be influenced by winter conditions more than by summer temperature (Harsch et al. 2009). Exceptions to these outcomes exist, however; for instance, where diffuse treelines have not shifted for 4,000 years (Millar et al. 2019) and where krummholz treeline ecotones did shift in response to summer warming during the late Holocene (LaMarche and Mooney 1972; Payette et al. 1989). Despite advances in understanding, the mechanisms that regulate specific treelines—particularly the influence of krummholz zones—relative to changing climates remain poorly understood (Körner 2012b; Camarero et al. 2017).

Stem growth in upper subalpine forest trees influences the response of alpine ecotones to changing climates. Controls on growth are complex and vary by species and location (Cook and Kairukstis 1990). For subalpine forests of semi-arid regions such as the American Great Basin, temperature has long been shown to control radial growth near upper treeline (UTL), whereas precipitation determined growth at the lower forest border (LaMarche 1974a, 1974b; LaMarche and Stockton 1974; Hughes and Funkhouser 2003). Recent studies, however, have revealed changing interactions of temperature and moisture on growth in trees near treeline that relate to contemporary climates (Salzer et al. 2009, 2014; Tran et al. 2017). In that zone, growth of bristlecone pines (Pinus longaeva) was related to micro-site, with precipitation playing a factor on warmer sites and temperature on cooler sites (Tran et al. 2017; Bunn et al. 2018). Moisture sensitivity appears to be a new factor affecting growth and potentially heralds changes in response at treeline in other species as well.

How this new pattern might relate to growth above treeline in krummholz has been little addressed. Studies of krummholz have focused on stem layering and genet age (King and Graumlich 1998), seedling recruitment (Daniels and Veblen 2004; Hofgaard, Dalen, and Hytteborn 2009; Schwab et al. 2016), and environmental controls of crown form, in particular the role of winter snowpack, wind abrasion, and low temperatures in flag release, survival, and mortality (e.g., Hadley and Smith 1987; Weisberg and Baker 1995; Smith et al. 2003; Millar et al. 2004; Devi et al. 2008; Maher, Nelson, and Larson 2020). Improved understanding of the influence of climate on stem growth in krummholz forms and its relation to upper subalpine forest growth could contribute to systematic understanding of treeline ecotone dynamics.

We sought to address specific uncertainties in response of treeline and in the krummholz zone with an integrated study of temperature and growth relationships in a representative North American treeline species, whitebark pine (WBP; Pinus albicaulis). WBP extends widely across western North American mountain ranges (Critchfield and Little 1966), where it is an important member of subalpine forest communities and commonly forms UTL. WBP has ready capacity to adopt stunted forms, but the expression of krummholz on the landscape varies geographically. In most parts of the species range the zone is narrow, occurring as a fringe above UTL or with scattered krummholz islands above treeline (Arno 1989). Throughout its range, WBP is a keystone ecological species, supporting critical ecosystem functions and serving as important wildlife habitat (Tomback, Arno, and Keane 2001). In many parts of its range WBP’s persistence is threatened by nonnative white pine blister rust (Cronartium ribicola), native bark beetles (Dendroctonous ponderosae), and effects of fire suppression (Keane et al. 1990). The impacts of combined disturbances on WBP mortality have been severe enough that the species has been designated as a candidate for listing under the Endangered Species Act in the United States and has been listed as endangered in Canada.

WBP is the dominant conifer of subalpine zones in California. We studied WBP treeline and krummholz vegetation in the eastern Sierra Nevada, California, where the species regularly develops broad zones of krummholz. Although little studied anywhere in its range, the WBP krummholz zone can harbor individuals of great age relative to upright trees (King and Graumlich 1998) and high genetic diversity (Rogers, Millar, and Westfall 1999); provide important wildlife habitat (Millar and Rundel 2016), snow retention, and watershed protection (Daly 1984); and serve as refugium in its apparent resistance to mountain pine beetle mortality (Maher 2018). Our study addressed the following questions:

1. What thermal conditions define the WBP forest-to-alpine ecotone in California, from treeline to upper krummholz (species) line (UKL)?

2. Does a consistent thermal threshold characterize UKL, and how does this differ from treeline temperatures?

3. How do treeline GSTs in WBP forests adjacent to extensive krummholz vegetation compare to global averages?

4. What relationships of radial growth to climate occur in krummholz WBP?

5. How do these compare to growth and climate relations of WBP upright tree growth at treeline?

Methods

Study area and field methods

We selected four sites along the eastern escarpment of the central Sierra Nevada, California (Figure 1). These comprised two pairs of northern and southern aspect slopes, one at Tioga Peak/Tioga Crest and the other at Dunderberg Peak/Virginia Canyon. WBP dominated the upper subalpine forest and krummholz zone at all locations, although occasional lodgepole pines (Pinus contorta) co-occurred. In this region, krummholz hedges and matts extended across as much as 400 m altitude, from ~3,100 m at UTL to ~3,500 m, where it formed upper species line. Study transects at each site extended from UTL into the alpine zone, with elevations varying by site, aspect, and sampling type (Table 1, Figure 2). All sites were underlain by metamorphic substrates.

Table 1. Study sites, transects, and plot locations for tree-ring and temperature measurements. UTL and the tree core plots lower are upright trees; the remaining plots are krummholz

		Ground distance (m)
Site (abbreviation)	Aspect	UTL-UKL	UTL-ALP	Slope inclination (°), UTL-ALP	Tree core plots	Temperature plots	Elevation (m)
Dunderberg Peak (SPL)	North	632	934	17		ALP	3,423
					0	UKL	3,338
					100		3,300
					200	MKR	3,265
					300		3,229
					400		3,198
					500	UTL	3,166
					600		3,135
Virginia Canyon (VRG)	South	424	598	28		ALP	3,534
					0^a	UKL	3,419
					100/200		3,375
					300^a		3,344
					400	MKR	3,309
					500^a		3,276
					600		3,244
					700	UTL	3,207
					800		3,179
Tioga Peak (TIO)	North	222	378	23		ALP	3,581
					0	UKL	3,308
					100	MKR	3,273
					200	UTL	3,238
					300		3,203
Tioga Crest (GAR)	South	771	1,023	13		ALP	3,450
					0	UKL	3,415
					100		3,378
					200	MKR	3,343
					300		3,309
					400		3,277
					500	UTL	3,243
					600		3,211

Note. ^aNo useable cores retrieved.

Figure 1. Map of the study region in the eastern Sierra Nevada, California, with location of four transect sites. See Table 1 for site names

Figure 2. Transect design. Transects extended from UTL to ALP, with intermediate plots at MKR and UK (upper species line). Two pairs of transects were installed on adjacent northern and southern aspects for a total of four transects

Temperature measurements

To assess thermal conditions, four plots were established along each transect, situated at consistent ecological positions (Table 1, Figure 2). These were chosen by field examination of vegetation and stand structure and included: UTL (conditions per Körner [2012a] and as described earlier), middle of the krummholz zone (MKR), UKL (upper species line, conditions per Körner 2012a), and the alpine zone above the upper species line (ALP). The MKR plot was located midway between the elevations of UTL and UKL, which varied by site; the ALP plot was situated above UKL by the same elevation interval (half the number of meters between UTL and UKL). Ground lengths of transects (UTL to ALP) ranged from 378 to 1,023 m (Table 1).

At each plot, a single 3-m PVC pole (“thermopole”) was anchored upright such that it extended at least 2 m above the ground. The bases of all but the ALP pole were in shade of forest or krummholz canopies, with the krummholz bases in deepest shade. Six temperature dataloggers (Maxim iButtons, models DS1921G and DS1922L) were deployed at each plot (Figure 3). The positions for the dataloggers were chosen to investigate different questions. Three units were wired onto the thermopole, at 2.0 m, 1.0 m, and 0.5 m height—these monitored air temperatures at different heights within the forest or krummholz canopy or, at the ALP plots, open conditions, including standard free-air temperature (2.0 m). A fourth iButton was buried 10 cm below the ground surface near the base of the pole, which is a proxy for 2.0 m surface air temperature daily means at the expected range of treeline summer temperatures (Körner and Paulsen 2004). Two iButtons were deployed on the ground surface, one at the pole and another in an opening within approximately 5 m of the pole; these were intended as indicators of snow cover duration. During the first 2.5 years we used 17-cm white PVC funnels as radiation shields for the pole units, securing the iButtons inside 8-cm white PVC “T” tubes and these within funnels. Because funnels broke off during the 2017 severe winter, we subsequently used T tubes alone, a design we have used in prior studies, and the protection we used for the ground and buried units from the start.

Figure 3. Plot design. At each plot, six temperature loggers were deployed: three on a PVC pole, at 2.0 m, 1.0 m, and 0.5 m height; one buried 10 cm at the base of the pole; and two ground surface loggers, one at the base of the pole and one approximately 5 m distant

IButtons were programmed for 4-hour intervals, starting at noon on the day of deployment. Dataloggers were installed between 7 June and 7 July 2015. All units were exchanged at least once annually, attempting to maintain continuous measurement, and uninstalled in October 2019.

Growth measurements

To assess radial growth in trees near UTL and in the krummholz zone we used the same four transects as in the thermal studies, adding plots to obtain approximately 30 to 35 m elevation intervals (following Salzer et al. 2014). Because the elevation extent of each transect differed, so did the number of coring plots at each site (Table 1). We labeled the uppermost location sampled as the “0 plot,” which coincided with the UKL thermopole plot. Tree ring plots downhill were labeled 100, 200, etc., and adjusted slightly so that a plot fell near the MKR and UTL thermopole plots. A final coring plot was approximately 25 to 30 m below UTL and represented typical upper-subalpine WBP forest conditions. Total number of tree ring plots per transect varied from four to eight (Table 1).

In July 2016 plots were selected in the field by designating the plot center as the location where our Global Positioning System device first registered the predetermined elevation for the plot within the transect. Initially we attempted to select sample stems randomly starting with the nearest stem to the plot center in the north direction and moving clockwise around the center to reach the desired sample number (sixteen). We maintained at least 5 m between sample stems to avoid coring multiple stems of tree clusters. Beyond that we choose stems for coring without regard to diameter except that stems for coring were ≥10 cm at coring height. When a tree did not yield an adequate core, we moved to the next stem. Radial cores were retrieved from sample pines with standard increment borers, including the pith as possible. At the UTL plot and the plot below UTL, cores were taken from upright tree stems at approximately 0.5 m height from the ground. In the krummholz zone, cores were extracted as near to the base of the central stem as possible. Often it was impossible to identify or access the base and/or to core in the dense krummholz canopy. In those cases cores were taken either from near-basal branches or along the stem within 2 m of the krummholz base.

Laboratory and statistical analyses

Temperature evaluations

Temperature data were downloaded from iButtons after retrieval from the field. After final retrieval in 2019, data were concatenated chronologically by site, transect, plot, and position for the extent of the study. Occasional gaps occurred due to device failures and loss or displacement of units. When this occurred, values were imputed using training data from the entire period of record for datalogger positions at that plot, estimated by the JMP Bayesian imputation procedure (JMP 12.2; SAS 2015). Approximately 8 percent of the entire data set was imputed.

Temperature data were summarized for the following indicators: annual mean daily temperature (T_ann); annual minimum mean daily temperature (T_ann-min); annual maximum mean daily temperature (T_ann-max); mean GST; mean growing season length (GSL), mean temperature of the warmest month (T_warm); mean number of snow-free days (DSF). GSL was calculated following Körner and Paulsen (2004): For the 10 cm buried position, the first day of the growing season was the spring date when the unit first exceeded 3.2°C (coinciding with weekly mean free-air temperature of 0°C) and the last day was the autumn date when the buried unit first reached 3.2°C; GSL was the number of days between these dates. GST was the mean of daily temperatures during the growing season as defined. For the 2.0 m air temperature estimates, GST was estimated following Paulsen and Körner (2014), using 0.9°C as the threshold for growing season. DSF was calculated by water year (WY) from the number of winter days for which the daily range of dataloggers was greater than 1.0°C (i.e., not snow covered; Lundquist and Lott 2008). To assess consistent temperature differences between the 2.0 m pole and 10-cm buried positions, we applied the matched pairs test in JMP (SAS 2015) based on mean daily temperatures.

Radial growth and climate analyses

Air-dried increment cores from the field were prepared for ring-width measurement, dating, and analysis using standard tree ring techniques (Holmes, Adams, and Fritts 1986; Cook and Kairukstis 1990). Cross-dating and quality control assessments were done using COFECHA (Holmes, Adams, and Fritts 1986; Krusic 2012). The chronology was built first by using series of well-correlated mature upright trees from the upper subalpine forest sites and adding the more recalcitrant krummholz series as the chronology became robust with samples. We used all intact cores to estimate tree and krummholz ages; for series that did not cross-date to the master chronology, we estimated tree age by counting rings (all trees were alive at the time of collection). We thus had two data set categories, one that composed the master chronologies (for growth correlations and climate analysis) and one for demography. In developing the master chronologies, we included series that had correlations ≥0.4. For both data sets, we grouped the samples into three sets: treeline (Table 1; tree core plots at UTL and below), middle krummholz (plots at UTL to MKR), and upper krummholz (plots at MKR to UKL).

For master chronologies, we imported raw ring widths into ARSTAN 44hc (Cook and Krusic 2016) for detrending and standardization of ring widths by biweight robust averaging to remove exogenous factors such as density and to maximize correlations with climate (Cook and Kairukstis 1990). Although the RCS (regional curve standardization) detrending algorithm (Esper, Cook, and Schweingruber 2002) is generally not used for short series, we used the method because it revealed low-frequency trends, showed increasing growth trends, and improved synchrony among the krummholz and tree sets. Data were detrended by a spline curve at 66 percent n (the period of record) at 50 percent cutoff. We also stabilized variance at 90 percent n with a 50 percent cutoff. Chronologies were developed for the three elevation groups at each of the four transect locations.

For assessment of radial growth with twentieth to twenty-first century climate, instrumental climate records were compiled from four Historical Climate Network (HCN) weather stations that we have used in past studies in the region: Lake Tahoe, Yosemite, and Independence, California, and Mina, Nevada (NCDC 2020). Following the approach of Millar et al. (2012), we statistically combined data from the individual stations into composite records (1899–2018), taking the first principal component for mean monthly minimum and maximum temperature, respectively, and annual WY precipitation. To test relationships of climate and ecological responses, we analyzed simple linear correlations by cross-correlations in the time series platform in JMP (SAS 2015). For the latter, we conducted a second-order least squares response surface model (JMP; SAS 2015) of standardized ring widths versus minimum and maximum annual temperature, respectively, and WY precipitation. We evaluated the behavior of these variables in second-order response surface models of the form (x + y + …) + (x + y + …)², where redundant interaction terms were omitted. Interactions in the surface model tested heterogeneity of slopes between the effects. Lead year cross-correlations were assessed in the time series analysis platform in JMP (SAS 2015). Two-year lead effect minimum and maximum temperatures and WY precipitation were the final variables used in the model.

Results

Temperature

Despite gaps in the record, mostly filled by imputation, temperature data were retrieved from all sites and plots for the 4.5 years of study, extending from early summer 2015 to autumn 2019 (Figure 4, Table 2). At UTL mean annual temperature was 3.0°C for 2.0 m (standard air temperature) and 3.6°C for 10 cm buried, with intermediate values at the 1.0 m and 0.5 m positions. Treeline growing season values were estimated from loggers at UTL plots at 2.0 m and 10 cm buried positions. These indicated mean GST 9.4°C across years and transects at 2.0 m, with mean GSL of 148 days and mean GST of 9.2°C with GSL of 128 days for the buried position. The longest snowfree duration occurred at 2.0 m (mean, 274 days) and shortest at the buried position (mean, 198 days; flat lines in Figure 4). High WY precipitation in 2017 resulted in lowest DSF for many transects, but despite variable precipitation in other years, DSF was not greatly affected.

Table 2. Temperature and growing season means based on daily means. Annual temperature (T_ann); minimum (T_ann-min) and maximum (T_ann-max) annual temperature; GST; GSL; warmest month temperature (T_warm); number of DSF

				Temperature (°C)
Site	Plot	Aspect	Position	Tann	SD	Tann-min	Tann-max	GST (°C)	GST SD	GSL (days)	Twarm (°C)	DSF (days)	DSF SD
SPL	UTL	N	Pole 2.0 m	2.4	7.4	−16	24	9.2	4.7	131	12.9	234	53
VRG	UTL	S	Pole 2.0 m	4.5	8.2	−18	21	9.3	5.1	172	14.3	364	1
TIO	UTL	N	Pole 2.0 m	2.3	7.5	−17	18	9.3	4.4	142	12.7	251	44
GAR	UTL	S	Pole 2.0 m	2.9	7.3	−17	23	9.6	4.5	146	12.6	248	59
Mean				3.0	7.6	−17	22	9.4	4.7	148	13.1	274	39
SPL	UTL	N	Pole 1.0 m	2.7	6.8	−17	24					224	49
VRG	UTL	S	Pole 1.0 m	3.8	8.3	−18	19					364	1
TIO	UTL	N	Pole 1.0 m	2.7	6.7	−22	19					221	49
GAR	UTL	S	Pole 1.0 m	3.5	6.2	−16	21					186	31
Mean				3.2	7.0	−18	21					249	32
SPL	UTL	N	Pole 0.5 m	2.5	7.1	−17	26					210	45
VRG	UTL	S	Pole 0.5 m	4.0	8.3	−18	20					365	1
TIO	UTL	N	Pole 0.5 m	2.4	6.7	−17	19					200	59
GAR	UTL	S	Pole 0.5 m	3.8	5.9	−17	19					169	32
Mean				3.2	7.0	−17	21					236	34
SPL	UTL	N	Buried 10 cm	3.7	5.6	−5	18	11.2	3.7	103	13.6	128	14
VRG	UTL	S	Buried 10 cm	3.7	5.8	−11	16	8.5	3.0	156	11.1	352	15
TIO	UTL	N	Buried 10 cm	2.9	5.0	−7	26	9.0	3.3	120	10.9	165	29
GAR	UTL	S	Buried 10 cm	3.9	3.8	−1	13	8.1	2.3	133	9.9	147	26
Mean				3.6	5.1	−6	18	9.2	3.1	128	11.4	198	21
SPL	UTL	N	Ground surface at pole	3.7	6.2	−6	20					146	15
SPL	UTL	N	Ground surface at distance	3.7	5.9	−4	21					124	39
VRG	UTL	S	Ground surface at pole	4.5	7.8	−13	20					364	2
VRG	UTL	S	Ground surface at distance	7.7	7.5	−7	23					297	31
TIO	UTL	N	Ground surface at pole	3.2	5.6	−7	17					141	15
TIO	UTL	N	Ground surface at distance	3.5	6.5	−8	20					128	22
GAR	UTL	S	Ground surface at pole	3.8	5.1	−4	17					150	20
GAR	UTL	S	Ground surface at distance	5.7	7.2	−4	23					168	29
Mean				4.5	6.5	−7	20					190	22
SPL	MKR	N	Pole 2.0 m	3.1	7.7	−18	19	9.7	4.6	156	13.5	248	84
VRG	MKR	S	Pole 2.0 m	4.2	8.0	−18	19	9.0	4.9	170	13.7	365	1
TIO	MKR	N	Pole 2.0 m	2.5	7.5	−17	20	9.3	4.6	147	12.6	251	95
GAR	MKR	S	Pole 2.0 m	2.8	7.8	−20	18	8.8	4.6	154	12.7	314	102
Mean				3.2	7.8	−18	19	9.2	4.7	157	13.1	295	71
SPL	MKR	N	Pole 1.0 m	3.3	7.9	−18	22					200	30
VRG	MKR	S	Pole 1.0 m	4.2	8.3	−18	20					365	1
TIO	MKR	N	Pole 1.0 m	2.2	7.3	−17	20					236	105
GAR	MKR	S	Pole 1.0 m	3.1	8.1	−19	20					307	74
Mean				3.2	7.9	−18	21					277	53
SPL	MKR	N	Pole 0.5 m	3.0	6.6	−11	20					157	25
VRG	MKR	S	Pole 0.5 m	3.2	8.3	−18	19					363	1
TIO	MKR	N	Pole 0.5 m	2.9	6.3	−10	20					186	58
GAR	MKR	S	Pole 0.5 m	3.9	6.8	−8	20					278	88
Mean				3.3	7.0	−12	20					246	43
SPL	MKR	N	Buried 10 cm	3.6	5.0	−2	16	10.2	3.0	116	12.1	139	26
VRG	MKR	S	Buried 10 cm	3.8	5.0	−8	16	8.8	3.0	142	11.2	194	19
TIO	MKR	N	Buried 10 cm	2.3	3.9	−4	12	7.6	2.2	103	8.4	173	115
GAR	MKR	S	Buried 10 cm	3.8	4.7	−4	16	9.1	2.9	128	11.2	157	30
Mean				3.4	4.7	−5	15	8.9	2.8	122	10.7	166	48
SPL	MKR	N	Ground surface at pole	4.7	6.2	−5	20					146	17
SPL	MKR	N	Ground surface at distance	4.3	7.1	−7	23					152	21
VRG	MKR	S	Ground surface at pole	3.6	6.6	−12	20					256	39
VRG	MKR	S	Ground surface at distance	n/a	n/a	n/a	n/a					n/a	n/a
TIO	MKR	N	Ground surface at pole	3.3	6.3	−8	20					148	34
TIO	MKR	N	Ground surface at distance	4.1	7.0	−10	21					144	24
GAR	MKR	S	Ground surface at pole	4.1	7.0	−9	23					174	35
GAR	MKR	S	Ground surface at distance	5.2	6.4	−4	21					154	22
Mean				4.2	6.7	−8	21					168	27
SPL	UKL	N	Pole 2.0 m	2.4	8.4	−19	18	8.8	4.6	158	12.7	356	18
VRG	UKL	S	Pole 2.0 m	2.4	7.6	−19	18	7.8	4.5	216	12.6	364	1
TIO	UKL	N	Pole 2.0 m	2.2	8.0	−20	19	9.1	4.6	149	12.8	297	85
GAR	UKL	S	Pole 2.0 m	1.6	7.3	−21	18	8.8	4.0	147	11.6	212	36
Mean				2.2	7.8	−20	18	8.6	4.4	168	12.4	307	35
SPL	UKL	N	Pole 1.0 m	2.3	8.3	−19	23					327	57
VRG	UKL	S	Pole 1.0 m	2.7	7.9	−20	19					365	0
TIO	UKL	N	Pole 1.0 m	2.1	7.6	−17	18					247	80
GAR	UKL	S	Pole 1.0 m	3.0	6.6	−14	21					159	17
Mean				2.5	7.6	−18	20					275	39
SPL	UKL	N	Pole 0.5 m	2.5	7.2	−11	23					206	41
VRG	UKL	S	Pole 0.5 m	3.8	6.5	−13	18					316	29
TIO	UKL	N	Pole 0.5 m	2.9	7.0	−14	19					200	69
GAR	UKL	S	Pole 0.5 m	3.3	6.3	−14	20					155	17
Mean				3.1	6.8	−13	20					219	39
SPL	UKL	N	Buried 10 cm	3.0	4.9	−3	14	9.3	2.7	112	11.0	129	16
VRG	UKL	S	Buried 10 cm	4.3	5.0	−3	16	9.1	3.2	151	11.6	226	59
TIO	UKL	N	Buried 10 cm	2.6	4.1	−4	13	8.1	2.2	106	9.2	152	76
GAR	UKL	S	Buried 10 cm	4.1	4.9	−4	18	9.8	3.4	124	11.6	148	21
Mean				3.5	4.7	−4	15	9.1	2.9	123	10.9	164	43
SPL	UKL	N	Ground surface at pole	3.6	6.0	−4	18					149	33
SPL	UKL	N	Ground surface at distance	4.4	7.1	−9	23					162	63
VRG	UKL	S	Ground surface at pole	4.5	6.0	−7	18					269	68
VRG	UKL	S	Ground surface at distance	4.8	7.6	−13	20					333	14
TIO	UKL	N	Ground surface at pole	2.8	5.5	−8	16					170	49
TIO	UKL	N	Ground surface at distance	3.2	6.4	−9	22					161	47
GAR	UKL	S	Ground surface at pole	4.8	6.3	−7	21					153	21
GAR	UKL	S	Ground surface at distance	4.9	7.1	−5	24					146	17
Mean				4.1	6.5	−8	20					193	39
SPL	ALP	N	Pole 2.0 m	1.5	7.7	−17	17	8.6	4.6	144	12.3	n/a	n/a
VRG	ALP	S	Pole 2.0 m	1.8	7.6	−21	18	8.2	4.1	149	11.5	316	39
TIO	ALP	N	Pole 2.0 m	1.6	7.9	−20	17	8.6	4.3	153	12.2	322	55
GAR	ALP	S	Pole 2.0 m	1.6	7.8	−18	16	8.6	4.1	150	11.6	308	75
Mean				1.6	7.8	−19	17	8.5	4.3	149	11.9	315	56
SPL	ALP	N	Pole 1.0 m	1.5	7.8	−20	18					282	18
VRG	ALP	S	Pole 1.0 m	2.0	8.1	−20	18					356	16
TIO	ALP	N	Pole 1.0 m	1.7	7.9	−19	17					293	78
GAR	ALP	S	Pole 1.0 m	2.1	6.9	−17	16					211	22
Mean				1.8	7.7	−19	17					286	34
SPL	ALP	N	Pole 0.5 m	2.6	6.8	−19	18					214	42
VRG	ALP	S	Pole 0.5 m	2.7	8.5	−22	19					348	33
TIO	ALP	N	Pole 0.5 m	1.7	7.9	−19	19					263	73
GAR	ALP	S	Pole 0.5 m	2.8	6.1	−13	16					186	22
Mean				2.5	7.3	−18	18					253	43
SPL	ALP	N	Buried 10 cm	2.9	6.0	−8	19	10.6	3.1	108	12.2	155	36
VRG	ALP	S	Buried 10 cm	4.7	7.1	−9	20	11.3	4.3	154	15.0	276	90
TIO	ALP	N	Buried 10 cm	2.2	6.5	−11	17	10.7	3.3	106	12.8	153	36
GAR	ALP	S	Buried 10 cm	5.0	5.9	−1	19	11.4	4.1	n/a	n/a	n/a	n/a
Mean				3.7	6.4	−7	19	11.0	3.7	123	13.3	195	54
SPL	ALP	N	Ground surface at pole	2.6	7.5	−14	24					184	58
SPL	ALP	N	Ground surface at distance	2.5	7.6	−15	22					225	78
VRG	ALP	S	Ground surface at pole	4.3	7.9	−11	21					316	7
VRG	ALP	S	Ground surface at distance	3.3	8.1	−15	29					348	28
TIO	ALP	N	Ground surface at pole	2.3	7.6	−14	20					180	51
TIO	ALP	N	Ground surface at distance	2.4	7.5	−14	20					192	59
GAR	ALP	S	Ground surface at pole	4.5	6.3	−5	19					151	21
GAR	ALP	S	Ground surface at distance	4.6	6.3	−6	20					149	11
Mean				3.3	7.4	−12	22					218	39

Note. n/a = insuffient data. See Table 1 for site names.

Figure 4. Overlay plots for mean 2.0 m air temperature and mean 10 cm buried temperature, July 2015 to October 2019 for paired TIO (north) and GAR (south) transect sites, upper treeline and upper species line. (a) UTL at TIO; (b) UKL at TIO; (c) UTL at GAR; (d) UKL at GAR

Comparing UTL conditions to UKL across years and transects indicated cooler mean annual temperatures at UKL, measured at 2.0 m as 2.2°C but not different at 10 cm buried (Table 2). Differences in air temperatures (0.8°C) and mean elevation (283 m) between UTL and UKL give a lapse rate of −2.8°C/1,000 m. The UKL loggers indicated cooler growing season thresholds at UKL than UTL. Mean GST estimated at 2.0 m was 8.6°C (difference of 0.8°C from UTL) with GSL of 168 days; at 10 cm buried GST was 9.1°C (difference of 0.1°C) with GSL of 123 days. At UKL, DSF for 2.0 m was longer than UTL (307 days), whereas for 10 cm buried DSF was shorter at UKL (164 days) than at UTL.

Above treeline, the ALP plots indicated strongest differences between buried and air temperatures of all elevations (Table 2). Mean annual temperatures were colder than other elevations for 2.0 m air temperature (1.6°C), whereas the buried position mean was warmest (3.7°C). The 2.0 m position had cooler mean GST (8.5°C) than other elevations, although the GSL (149 days) was not different from UTL. By contrast, mean GST for the 10 cm buried position was much warmer than other elevations (11.0°C), with shorter GSL (123 days). Snow cover differed between positions, with DSF estimated as 315 for 2.0 m vs. 195 days for the buried position.

Annual minimum mean daily air temperatures over years and transects at the 2.0 m position generally followed expected elevation patterns along transects, from warmest (low) to coolest (high): −17°C (UTL), −18°C (MKR), −20°C (UKL), and −19°C (ALP; Table 2). Corresponding values for annual maximum mean daily temperatures similarly followed the elevation gradient: 22°C (UTL), 19°C (MKR), 18°C (UKL), and 17°C (ALP). The warmest month mean temperatures for the 2.0 m position were 13.1°C (UTL), 13.1°C (MKR), 12.7°C (UKL), and 11.9°C (ALP). These differences translate to mean lapse rates of −1.2°C/1,000 m for warmest month, −8.3°C/1,000 m for T_ann-min, and −5.6°C/1,000 m for T_ann-max.

Aspect differences in temperature followed expectations in mean annual air temperature for the SPL (Silverpine Lake)/VRG (Virginia Canyon) pair, with northern aspects colder than southern in all cases (Table 2). For the other pair (TIO (Tioga Peak)/GAR(Tioga Crest)) and for other annual means, aspect differences were inconsistent. Snow persisted longer on northern aspects (156 DSF) than southern (224 DSF), creating mostly consistent differences between aspects. In WY 2017 and 2019, northern aspects had fewer DSF than other years, whereas southern aspects generally showed little relationship between DSF and WY. Overall, northern aspects had significantly lower DSF than southern aspects, with the VRG transect particularly dry (p < .05). For the important GST and GSL estimates, differences were apparent between aspects but inconsistent in some variables (Table 2). At UTL, differences in GST were small (range = 0°C–0.4°C) at 2.0 m with no aspect effect but larger for the buried estimate (range = 0.5°C–3.1°C), with northern aspects unexpectedly warmer than southern. Mean GSLs showed large differences, with northern aspects shorter (GSL 137 days) and southern aspects longer (159 days).

Inconsistent aspect patterns occurred at UKL (Table 2). Annual mean temperatures between aspects were identical for one pair of sites; for the other (TIO/GAR), the southern aspect was 0.6°C cooler. Mean daily minimum and maximum annual temperatures did not show an aspect effect. Threshold GSTs measured at 2.0 m were warmer on the northern aspects by 0.3°C to 1.0°C relative to their paired southern sites. GSL was shorter by fifty-eight days on one northern aspect relative to the southern but basically the same number of days for the other north–south pair. At the 10-cm buried positions, GST was warmer by 1.7°C on the northern aspect than on the southern aspect in one pair of transects and cooler by 0.2°C in the other. GSL at the buried positions was shorter on northern aspects than on southern aspects by twelve to forty-five days.

Radial growth

A total of 456 increment cores composed the demography data set (Figure 5, Table 3). Overall the trees were relatively young, and series ranged in age (estimated from ring number) from 28 to 475 years, with a mean of 146 years. Pooled by elevation and crown form strata, stems in the upper elevation krummholz zone were youngest (mean 99.6 years), and stems in the middle elevation krummholz and the lower elevation forest zones were older and similar in mean age (158.2 and 166.5 years, respectively) although the forest zone had more older trees than the middle krummholz.

Table 3. Mean number of rings (estimate of age) and number of series (cores) contributing to tree-ring demography database for krummholz and tree plots

					Number of rings
Site	Aspect	Elevation zone	Stem form	Plots included	Mean	SD	Number of series
SPL	N	Upper	Krummolz	0, 100	124.0	43.3	40
VRG	S	Upper	Krummolz	100/200	86.0	18.9	15
TIO	N	Upper	Krummolz	0	94.2	30.5	16
GAR	S	Upper	Krummolz	0, 100	94.3	30.9	42
Overall					99.6	30.9	113
SPL	N	Middle	Krummolz	200, 300, 400	196.3	82.0	70
VRG	S	Middle	Krummolz	400, 600	146.4	49.3	52
TIO	N	Middle	Krummolz	100	173.6	60.6	15
GAR	S	Middle	Krummolz	200, 300, 400	116.5	43.5	56
Overall					158.2	58.9	193
SPL	N	Lower	Tree	500, 600	225.6	58.1	44
VRG	S	Lower	Tree	700, 800	144.6	76.0	31
TIO	N	Lower	Tree	200, 300	167.3	72.7	38
GAR	S	Lower	Tree	500, 600	128.6	31.0	37
Overall					166.5	59.6	150

Figure 5. Tree ring series, demography data set. Dated series from 456 live krummholz and upright trees sampled along four transects. Series pooled and arranged by upper elevation krummholz, middle elevation krummholz, and lower elevation (upper treeline) forest strata. All stems were live at the time of sampling

Ring-width behavior differed between the set of krummholz stems and that of forest trees. In general, trees had ring characteristics that are typical of WBP in Sierra Nevada environments (King and Graumlich 1998; Millar et al. 2012). That is, the ring series were relatively complacent, and generally fast growing with occasional periods of narrow rings. Filtering the master chronology to series with correlation ≥0.4 resulted in retaining 50 percent of the total number of tree series. The inability of 50 percent of the series to correlate at higher value was somewhat surprising, in that the cores were good quality and series excluded from the master chronology were not removed due to problems such as rot, knots, or poor cores. Nonetheless, based on our prior work with WBP dendrochronology, this was not unusual.

What was unusual in our experience was the far lower correlation of krummholz series to each other and to the master chronology, as estimated in COFECHA. These cores also were high quality, mostly fast growing, and free of damage that would prevent accurate dating. In some cases, series had intervals of extremely narrow rings; mostly these were asynchronous across the krummholz set (as frost rings might be), although neighboring stems within a plot occasionally exhibited coincident behavior. As a result, far fewer krummholz than tree series correlated with its master chronology, with only 19 percent of the total krummholz series correlated ≥0.4. No obvious differences in aspect characterized the ability to correlate the series, nor did age affect dating of one set versus the other.

The master chronology comprised 136 series. Detrending by the RCS method indicated synchronies in growth over the last 120 years among the krummholz and tree groups (Figures 6a–6c). Large declines in growth occurred during 1920–1935; growth was greatest in the 1940s. The crown form groups showed simultaneous growth declines during 2000–2010, with increases in growth in the past decade, although the tree group showed a plateau and decline in growth in the last few years. A positive centennial trend in growth was present for both sets of krummholz stems but not for the trees.

Figure 6. (a)–(c) Whitebark pine radial growth and (d)–(f) composite HCN climate station means (monthly means, averaged annually, 1895–2020). Climate data averaged for Tahoe, Yosemite, Mina, and Independence HCN stations. Black dots, ring widths; red line, mean ring widths; green line, linear fit; blue curve, spline fit at lambda 10

Model correlations for climate variables in the regression analyses for ring width (Table 4A) indicate that trees responded significantly if weakly to maximum temperature (negative correlation, p < .05) and WY precipitation (positive, p < .05). By contrast, the krummholz groups (the small set that correlated with the master chronology) had different and stronger correlations with climate, responding positively to increases in minimum temperature (p < .0001). Over the past 120 years, maximum temperature and WY precipitation (Figures 6d and 6f) showed variability but no strong century trends. Increases in temperature in recent years correlated with declines in tree growth. The strong positive centennial increase in minimum temperature (Figure 6e) correlated with increased growth in krummholz. Aspect differences in correlations further distinguished krummholz from tree groups, with trees having the stronger aspect effect (Table 4B), although overall trees were more complacent with regard to climate (Table 4A). Interestingly, series from the uppermost north-facing plots had a stronger relationship with WY precipitation than those on south-facing slopes.

Table 4. Model correlations and climate variables included in the best-fit stepwise regression analyses for radial ring-width models, based on the date range 1895–2020.^a

		Trees	Middle krummholz		Upper krummholz
A. Crown form (elevation) groups
R^2		0.15	0.34		0.29
Model		0.0376	<0.0001		<0.0001
Tmax*		−, 0.0469
Tmin*			+, <0.0001		+, <0.0001
WY precipitation*		+, 0.0210
	Northern			Southern
	Trees	Middle krummholz	Upper krummholz	Trees	Middle krummholz	Upper krummholz
B. Aspect and crown form (elevation) groups
R^2	0.17	0.31	0.26	0.27	0.15	0.28
Model	0.02	<0.0001	<0.0001	<0.0001	0.0026	<0.0001
Tmax*	−, 0.28	−, 0.41	−, 0.0219	—	—	—
Tmin*	−, 0.24	+, <0.0001	+, <0.0001	+, <0.0001	+, 0.0007	+, <0.0001
WY precip*	+, 0.38	+, 0.0012	+, 0.0121	+, 0.011	+, 0.09

Notes. ^a(A) All samples, grouped by crown form (elevation). (B) Samples grouped by aspect and crown form. Climate associations are derived only from tree and krummholz series that had correlations >0.4 to the master chronology.

*Sign of the correlation; probablility.

Discussion

Potential biases in data

We report the first investigation of thermal relationships of WBP at treeline and within the WBP krummholz zone and, to our knowledge, the first study of threshold GSTs and lengths in krummholz zones for any conifer species. In that we sought to evaluate differences for krummholz vs. forest trees within our study area and to compare with treeline results from studies elsewhere, it is important that we review potential biases. Temperature is difficult to measure accurately under field conditions, in particular to protect loggers from direct sun effects. Our air units (loggers on poles), although protected by white PVC containers, potentially experienced heating in some cases. We summarize our assessments of bias here; more detail is provided in the Supplemental Material.

In the krummholz plots, the 0.5 m and 1.0 m pole units were consistently in deep shade of dense crowns and likely recorded temperatures relatively accurately. All pole positions at UTL plots similarly had dense overstory shade, whereas the krummholz 2.0 m pole positions mostly extended above the krummholz crown and were exposed to direct sunlight. Buried units in all plots except ALP were under deep shade and thus presumably were recording standard air temperatures when snow free with relative accuracy (Körner and Paulsen 2004). To assess effects of direct solar exposure on the upper pole units, we compared June to September temperatures between 2.0 m and buried units at UTL and UKL transects and found the pole units to record mostly warmer values (Supplemental Material). The buried units, however, were variously covered by snow during these months, biasing their values toward cooler temperatures. To avoid problems of comparing relative values within our study, we used the approach and results of Strachan and Daly (2017) to provide an external assessment. In that comparison, our 2.0 m pole units gave no evidence for bias either warm or cold (Supplemental Material).

Thermal characteristics of upper treeline and species line; response to climate change

A striking finding in our 4.5-year study was the high GST at treeline (mean 9.3°C), which was 2.6°C higher than the global treeline threshold of 6.7°C calculated by Körner and Paulsen (2004; 10 cm buried estimate) and more than three times beyond their ±0.8°C standard deviation. Considering Körner and Paulsen’s (2004) finding that temperate and Mediterranean treelines have mean GST values higher (7°C–8°C) than other biomes, our GST estimates were warmer even than these. The two California sites (Tioga Pass, Sierra Nevada) used in Körner and Paulsen’s (2004) study, which lie only 7 km from our nearest transect, yielded means warmer than the other Mediterranean locations in their global study. The authors did not include these values in estimating the global means due to what they felt was a bias effect from snowpack, explained below. GST estimates from their Tioga Pass sites were 9.4°C (southern aspect) and 8.0°C (northern aspect), which closely resemble our estimates. Our GSL (143 days), estimated using Körner and Paulsen’s (2004) buried logger method, was longer than their global threshold of 94 days but similar to their global mean GSL of 135 days. A long growing season is expected relative to other parts of the world, given the extended Mediterranean-climate summers in the Sierra Nevada.

Körner and Paulsen (2004) interpreted their exceptionally warm GST and short GSL at Tioga Pass to result from a single anomalous year of measurement. That year (WY 1997) included a very wet winter, with snowpack lasting into summer (DWR 2020). When buried dataloggers remain under snow (and thus record temperatures below the growth threshold) in summer, the estimated GST is elevated and GSL biased shorter. The five summers of our study, however, followed winters that spanned the range of California snowpack depths, from very heavy years (WYs 2017 and 2019) to exceptionally low (WY 2015) and low–moderate (WYs 2016 and 2018) years. Pooling data across years as we did should not only avoid effects of a single extreme year but provide more biologically meaningful values than data from any single year.

We had a second method, namely, air temperature, for estimating GST and comparing with the buried loggers. Temperature measured at 10 cm below the soil surface at our latitude is a surrogate for mean daily air temperature (2.0 m) during the growing season (Körner and Paulsen 2004). Following Paulsen and Körner’s (2014) methods for air temperature assessment, our GST was again warmer, in this case by 2.9°C, than the global mean for air-estimated GST of 6.4°C (Körner 2012a; Paulsen and Körner 2014). Combining evidence, we propose that the treeline GST and GSL values estimated in the present study reflect relatively accurate values for the Sierra Nevada, suggest that those in Körner and Paulsen (2004) were also not significantly biased, and argue that Sierra Nevada WBP treeline threshold values differ from other situations. We take up potential reasons for these deviations after further consideration of treeline and krummholz conditions.

To our knowledge there are only two published accounts of treeline growing season thresholds for pines in the Great Basin. Bruening et al. (2017) reported thresholds based on air temperature for Great Basin bristlecone pine and foxtail pine (Pinus balfouriana). Their field methods were similar to ours, although they used loggers hung at 1.0 m in trees and measured temperatures only in one year (2015). Further, before calculating temperature thresholds, including GST (J. Bruening, pers. comm. 21 February 2020), 1.5°C was subtracted from all raw temperature measurements to adjust for disequilibrium of present-day conditions from Little Ice Age (LIA) effects. For Great Basin bristlecone pine in the Snake Range, eastern Nevada, and White Mountains, California, Bruening et al. (2017) reported GST as 7.2°C and 6.3°C, respectively, and for foxtail pine in the southern Sierra Nevada, California, they reported GST as 5.5°C. Realigning their values (adding 1.5°C) to compare with the present study gives respective GSTs of 8.7°C, 7.8°C, and 7.0°C. Though these values are cooler than our estimates for WBP in the Sierra Nevada, they are warmer than the global mean threshold given by Paulsen and Körner (2014). In the Bruening et al. (2017) study, GSL ranged from 147 to 153 days, which is similar to our estimates. The other study of treeline temperatures in the Great Basin investigated a southern aspect limber pine (Pinus flexilis) ecotone in the arid Toquima Range, central Nevada (Millar in press). Treeline attained at 3,311 m and GST estimated from air temperatures (methods as in the present work) across nine years (2009–2018) was 12.2°C with GSL of 106 days. This is 2.8°C warmer than the mean GST for WBP.

Thermal relationships we report for the krummholz zone do not have relevant field comparisons in the published literature but differ from theoretic expectations. The low, shrub-like stature of krummholz is expected to experience different microclimates from conditions of subalpine forest trees, whose crowns are coupled to atmospheric air temperatures (Körner et al. 2016). Decoupling in the krummholz vegetation should result in crown and ground surface temperatures warmer than under forest crowns at treeline (Hadley and Smith 1987; Grace 1988; Körner 2007). In our study, however, mean annual temperatures and GST followed patterns expected from atmospheric lapse processes; that is, cooling with increased elevation. We interpret the dense krummholz crowns to provide consistent and effective shading, which resists solar radiation–induced microclimatic variation on the ground and air surface below the crowns. Further, our 2.0 m air temperature measurements, which showed expected elevation-related decreases, do not represent foliar temperature, which likely would be much different (Grace 1988). Notably, the GST at upper species line (8.8°C, 0.5°C cooler on average than treeline GST) is also considerably warmer than the range of values given for global treelines (Körner and Paulsen 2004), suggesting, at first blush, that current thermal conditions at upper species line in the Sierra Nevada are appropriate to sustain upright tree growth.

We turn to consider potential reasons for these anomalous results; that is, GST warmer than expected at treeline, cooler temperatures (rather than warmer) through the krummholz zone to upper species line, and GST at upper species line apparently warm enough to support trees. Our results are similar to other Sierra Nevada and Great Basin studies that found warm treeline GSTs, as described above. Further, the finding of very warm upper species line GST (11.1°C) for limber pine in the Toquima Range (Millar in press), also higher than global treeline values, corroborates the range of values we estimated. One possibility for these findings follows explanations given by Körner and Paulsen (2004): In general they found that Mediterranean-climate treelines form at consistently warmer GSTs than do those in other biomes. Some treelines in Mediterranean regions, however, were not determined by climate but by taxon-specific factors. In that they found WBP treeline in the Sierra Nevada to occur at an isotherm much warmer than expected (yet similar to our results), unique characteristics related to WBP in the Sierra Nevada might limit the species ability to ascend to higher elevations. In particular, we propose one possibility that WBP treeline in the Sierra Nevada is influenced as much or more by growing season moisture availability than temperature. The limiting role of soil moisture, exacerbated under warming climates, has been demonstrated previously for Great Basin subalpine species in seedling and tree establishment (Millar et al. 2015, 2018, 2019). Further, the positive response of radial growth to increased moisture found in the current study as well as prior Great Basin pine studies (Millar et al. 2004, 2012, 2018, 2019; Tran et al. 2017; Bunn et al. 2018) points to a potential role of precipitation limitation for forests in these long, dry Mediterranean summer climate regions. The role of growing season moisture restrictions was considered but dismissed by Körner and Paulsen (2004) and deserves further attention in these warm, semi-arid treeline regions.

Beyond these considerations, WBP treelines in the Sierra Nevada might be in disequilibrium with current climates and lagging in response to current temperature. The global LIA brought persistent cold and wet conditions to the Great Basin from 1400 CE to 1920 CE, initiating the largest glacial advance in the Sierra Nevada since the late Pleistocene (references in Millar et al. 2019). Decreases in summer minimum temperature of 0.2°C to 2°C and increases in winter precipitation of 3 to 26 cm compared to modern times have been interpreted to explain glacial conditions (Bowerman and Clark 2011). Summers significantly cooler than present persisted to 1930 CE (Feng and Epstein 1994). Five centuries of cold conditions could have depressed tree and species lines or at least enforced severe krummholz conditions. Depression of our GST values by 1.5°C to 2°C would put LIA treeline GST within the range expected for Mediterranean conditions elsewhere in the world. As noted, Bruening et al. (2017) accepted an LIA disequilibrium in their analyses for bristlecone and foxtail pine treelines in the Great Basin, adjusting their observations accordingly. Other studies have interpreted anomalous treeline demography and temperatures to result from lag with contemporary climate and especially disequilibrium as a result of LIA conditions (Payette et al. 1989; Rupp, Chapin, and Starfield 2001; Harsch et al. 2009; Körner 2012a; Alexander et al. 2018; Smithers et al. 2018; Malanson et al. 2019).

The possible role of disequilibrium is informed by the age distribution of stems in our plots. Approximately 25 percent of stems in subalpine forest and middle krummholz plots and approximately 75 percent of stems in the highest krummholz plots, established since the end of the LIA (approximately 1920–1930 CE; Feng and Epstein 1994; Bowerman and Clark 2011), when summer temperatures began to warm consistently. Prior studies have reported the response of Sierran conifers to warming of the past 100 years to be densification of existing subalpine and krummholz communities without treeline advance (Vale and Vale 1994; Millar et al. 2004; Dolanc, Thorne, and Safford 2013). Although large-diameter stems grew in our krummholz zones, their occurrence was rare enough that by chance very old individuals were not included in our plots. By contrast, in a study of long-lived krummholz WBP near our study transects, King and Graumlich (1998) targeted large-diameter stems for dating. They reported stems extending throughout the elevation range of the current krummholz zone as more than 1,000 years old and likely as old as 2,000 years. The occurrence of such individuals points to the ability of krummholz WBP to persist through great climate variability, including not only the cool LIA centuries but the warm, dry centuries of the Medieval Climate Anomaly (Stine 1994; Meko et al. 2001). The large proportion of young stems at our highest krummholz plots corroborates a hypothesis of infill around preexisting old stems with warming of the past century. The young plants at high elevations could originate not only from seeds of existing old-aged krummholz individuals (Maher 2018) but from seeds brought by birds from lower subalpine forest or middle-elevation krummholz stands (Rogers, Millar, and Westfall 1996). Because we did not study seedling age or abundance within forest or krummholz plots or above upper species line, we cannot conclude whether recruitment is still ongoing at those elevations. Anecdotal observations, however, as well as prior studies (Vale and Vale 1994; Dolanc, Thorne, and Safford 2013) have revealed little evidence for seedling recruitment above UKL into the current alpine zone.

A final consideration is that the krummholz WBP zone in the Sierra Nevada acts as a buffer to changing climates, with increases in stem densities during warm climate periods and decreases during cool periods. This follows observations of other krummholz-forming conifers, whose treelines appear relatively unresponsive to climate change (Harsch et al. 2009), and points to the influence of tree form on treeline dynamics (Harsch and Bader 2011). In addition to a role for growing season moisture limitation, winter conditions might enforce maintenance of the krummholz zone (Holtmeier and Broll 2010; Harsch and Bader 2011). Winter freezing damage, wind, abrasion, and snow breakage might regulate stem and branch dieback and seedling establishment, as found in other studies with WBP (Daly 1984; Millar et al. 2004; Maher, Nelson, and Larson 2020). Despite relative persistence of krummholz zones, however, there is evidence that this vegetation community does react to climate change over time, at least at its upper and lower borders. LaMarche and Mooney (1972) observed a 100-m shift downward of the UKL during the neoglacial interval (4,000–2,000 years ago) in bristlecone pine in the eastern Great Basin, and Payette et al. (1989) documented shifts up in spruce krummholz in spruce (Picea) in Canada relating to historic climate. Thus, a buffering effect relative to changing plant density might be an early response to climate change.

Radial growth of tree form versus krummholz

Whereas thermal relations in the krummholz zone were similar to those of the treeline forest, with expected elevation-related differences, radial growth differed unexpectedly between the forms. Despite interannual weather variability, ring widths of mature trees near treeline were relatively complacent. Given that the sites were chosen for representation and not for unusual stress, this was not surprising (King and Graumlich 1998). A well-correlated chronology could be constructed from the tree series, which showed responsiveness to centennial climate. By contrast, the krummholz series, with very low correlations to each other or to the tree chronology, suggest that something other than regional climate influenced their growth.

One possibility is that elongating branches at the top of krummholz crowns experience and respond to strong microclimatic differences (Hadley and Smith 1987; Grace 1988; King and Graumlich 1998) as trees in different micro-sites near treeline do (Bunn et al. 2018). Although we measured air temperatures within and above the krummholz crowns, foliar temperatures are likely very different, with elevated day and depressed night extremes relative to free air (Grace 1988; Körner 2007), and these would affect radial growth. Another possibility is that the highly compact, dense krummholz crowns, with multiple rooting bases, create strong competition effects on growth of individual stems. If highly local competitive factors dominate over regional or even micro-site climate, then ring widths would not show high correlation to one another or to forest trees.

For the krummholz stems that could be included in the master chronology, their stronger and differing climate correlations from tree stems point to further differences in WBP crown-form growth. Strong correlations to minimum temperature in krummholz, and thus increases in growth as T_min increased over the past 120 years, conform to expectations about highest elevations being temperature sensitive. That expectation, however, has previously been evaluated for trees near treeline, not krummholz forms. In eastern Sierran WBP, the sensitivity to minimum temperature and precipitation resembles the situation of Great Basin bristlecone pine, where climate correlations near treeline varied by micro-site (Salzer et al. 2014; Bunn et al. 2018). This is corroborated by the aspect differences in WBP, wherein trees showed low climate sensitivities on northern aspects but strong correlations to minimum temperature and precipitation on southern aspects. By contrast, krummholz had strong correlations to both variables on northern and southern aspects.

Conclusions

A multiyear study of temperature across the extensive treeline–krummholz (shrub) ecotones of WBP in the central-eastern Sierra Nevada, California, indicated treeline GST of 9.3°C, which is 2.6°C warmer than mean global thresholds of Körner and Paulsen (2004) and Paulsen and Körner (2014). Declining temperatures with increasing elevation across the krummholz zone suggested that near-ground and free-air environments in krummholz were not experiencing expected warmer microclimates relative to forest conditions. At the UKL (upper species line), GSTs were 2.4°C warmer than the average global treeline threshold of Paulsen and Körner (2014), suggesting that by thermal criteria these elevations could at present support tree growth. Lack of upright tree establishment above treeline and in the krummholz zone might result from taxon-specific idiosyncrasies of WBP in the Great Basin Mediterranean-climate region, including growing season moisture limitations and winter stresses that enforce krummholz persistence. Disequilibrium with LIA temperatures likely also creates lag in response of treeline to warming temperatures of the last century. Overall, we propose that the krummholz form influences the response of UTL, the density of vegetation in the krummholz zone, and likely movement of upper species line. As a result, the WBP treeline ecotone in the Sierra Nevada appears at present to be relatively stable.

Radial growth in WBP trees and krummholz showed differing responses to climate of the last 120 years. Low correlations of ring widths among krummholz within and across sites suggest that many individuals are responding to nonclimatic conditions and may be more influenced by intraspecies competition and microclimate. Ring-width patterns of upright trees, though having low sensitivity to climate, were in the range expected for WBP. Of the series that did correlate to a climate-sensitive chronology, krummholz growth responded positively to minimum annual temperature, which rose over the last 120 years, whereas tree growth responded negatively to maximum annual temperature and positively to precipitation, neither of which had dominant centennial trends. Though trees near treeline showed microsite differences (by aspect), paralleling recent studies with Great Basin bristlecone pine, krummholz growth was not responsive at this scale.

Results from the present study underscore the variable nature of treeline response to climate change and add to the examples of mountain conifers that may be exhibiting lag effects and have not shifted with contemporary warming. The WBP krummholz zone appears to affect dynamics of treeline ecotones, acting as buffers to treeline movement, increasing in stem density during warming, and responding to different growth cues than upright trees near treeline. If WBP in other locations and/or other krummholz-forming conifers behave in a similar manner, assessment of dynamics under changing climates can be made more accurate by taking into account climate response by these different conifer forms and effects of regional climates.

Acknowledgments

We thank Christian Körner (University of Basel) for discussions over many years and especially for thoughtful and detailed conversations during our data analysis period. Discussions with other colleagues have shaped our thoughts and provided important insights; in particular, we thank Andy Bunn (Western Washington University), Greg Pederson (USGS), Scotty Strachan (University of Nevada, Reno), and Connie Woodhouse (University of Arizona). Chrissy Howell (USDA Forest Service) reviewed the draft manuscript. We thank Lynn Resler (Associate Editor, AAAR) and two anonymous reviewers for constructive review comments.

Disclosure statement

No potential conflict of interest was reported by the authors.

Supplementary material

Supplemental data for this article can be accessed on the publisher’s website.

Additional information

Funding

Operational funds from the USDA Forest Service supported this study.