Sufficient conditions for rapid range expansion of a boreal conifer
White and black spruce are the dominant conifers at Arctic treelines and the boreal forest–tundra ecotone more generally in North America, with white spruce dominating on better drained sites. White spruce reaches its northwestern-most limit in Alaska, USA, at 68.1º N, 163.2º W. For comparison, the northeastern range extent of the species26 is Labrador, Canada, at 57.9º N, 62.5º W (ref. 12), giving an east–west range of >100º in longitude. Of the approximately 6,500-km-long northern boundary of white spruce in North America, 10–15% is located in Alaska’s Brooks Range, where white spruce is the dominant treeline tree.Study areaThe 1,000-km Brooks Range is a high-latitude mountain range dividing Arctic tundra from boreal forest in Alaska. The mountains and nearby lowlands are notable for their wilderness character, protected as a near-contiguous conservation area of >150,000 km2. In the east between the Arctic Ocean’s Beaufort Sea and the uppermost Yukon River basin, the range is cold and dry, reaching 2,736 m above sea level. The south slope of the eastern Brooks Range is included in Alaska’s Northeast Interior climate division, where precipitation is among the lowest in the state51. Descending to the Chukchi Sea in the west, the range is included in Alaska’s West Coast climate division, where precipitation is the highest in northern Alaska51.The Noatak and Kobuk rivers flow in their entirety above the Arctic Circle, draining the western Brooks Range. Both rivers empty into the Chukchi Sea near Kotzebue, Alaska (Fig. 1a). The Baird Mountains of the southwestern Brooks Range separate the Kobuk from the Noatak basin, and the De Long Mountains of the northwestern Brooks Range separate the Noatak from the river basins of the North Slope and from the Wulik basin, located northwest of the Noatak basin. The lower basins of the Noatak and Kobuk rivers are included in the West Coast climate division, with greater precipitation, warmer winters and cooler summers than in the Central Interior climate division and greater precipitation and warmer temperatures than in the North Slope climate division51. The upper basin of the 700-km Noatak River lies at the intersection of all three climate divisions, which warmed from 1949 to 2012; December–January precipitation increased from 1949 to 2012 in the West Coast climate division, as did North Slope winter precipitation from 1980 to 2012 (ref. 52).The Noatak River basin is entirely protected within federal conservation units. Its vegetation includes dwarf, low and tall shrub tundra communities that cover about 60% of the 33,000 km2 basin53. Tussock sedge tundra covers another 30%, and wetlands and barrens cover most of the remainder. The main valley and tributaries along the lowest 200 km of the Noatak River support stands of white spruce, typically associated with a deeper active layer or an absence of permafrost. The treelines bounding these forests have long been identified as the northwest range extent of white spruce26.The upper Noatak basin, a 500-km reach, is underlain by extensive continuous permafrost54. It has been considered empty of spruce since US Geological Survey (USGS) geologist Philip Smith explored the Kobuk, Alatna and Noatak rivers by canoe in 1911 (ref. 55). The adjacent Kobuk and Alatna river basins support boreal forests of black and white spruce, paper birch and aspen along much of their lengths. By surveying transects at and beyond hydrological divides separating the Noatak, Wulik, Kobuk and Alatna river basins, as well as further east in the Brooks Range (Fig. 1a), and informed by very high-resolution satellite scenes (Fig. 1b and Supplementary Figs. 1–13), we documented the locations of over 7,000 individual spruce colonists (Extended Data Fig. 1b–d and Supplementary Figs. 1–3). Overall, we traversed 22° of longitude (141–163° W) in the field, mostly along the treeline from Canada to the Chukchi Sea, locating dozens of populations of colonizing spruce (Fig. 1a) above alpine and beyond Arctic treelines (see ‘Regional extent of colonization’).The primary AOI (Fig. 1a) included the USGS Hydrological Unit Code (HUC) 10 watersheds Kaluich, Cutler, Amakomanak and Imelyak located in the HUC 8 Upper Noatak Subbasin. However, we also documented (longitude, latitude, distance from established treeline) fast-growing, healthy spruce well beyond established treelines within six additional western Arctic watersheds, each separated by over 30 km in the western Brooks Range and 80–200 km distant from the AOI. These populations are within the far upper reaches of the Noatak basin (Lucky Six Creek, 67.594° N, 154.858° W; Kugrak River, 67.428° N, 155.723° W; Ipnelivuk River, 67.552° N, 156.293° W; upper Wrench Creek, 68.251° N, 162.617° W); 25 km northwest of the nearest established treeline and outside the Noatak basin in the Wulik River valley (68.120° N, 163.219° W); and along the Chukchi Sea coast (67.041° N, 163.114° W). We also note that, in the central Brooks Range, humans have actively or inadvertently disseminated spruce seeds and juveniles on the North Slope, with individual white spruce germinating and surviving there for at least 20 years37,56.Patterns of expansionDigitizing spruce shadowsWe used cloud-free Maxar Digital Globe WorldView-1 and WorldView-2 satellite scenes (WV; https://evwhs.digitalglobe.com/myDigitalGlobe/login) of snow-covered landscapes from three missions in early spring 2018, a near-record year for snow depth in northwest Alaska (Fig. 1b, Extended Data Table 1 and Supplementary Figs. 1–13). Ground sample distances of 0.47–0.5 m, a root-mean-squared error of 3.91–3.94 m and off-nadir angles of 5–25º with low sun-elevation angles of 18–27º provided clear images from which to digitize the lengths of individual spruce shadows and identify their locations (Supplementary Information sections 1.2 and 1.3). One technician (S. Taylor), supervised in quality assurance and quality control (QAQC) by R.J.D., digitized 5,986 shadows (densities in Extended Data Fig. 1b, locations in Supplementary Fig. 1) on GEP using WV images as super-overlays. The technician identified all spruce shadows across the imported image tiles and then digitized them as line segments from base to shadow tip.The super-overlays degraded the imagery somewhat, making small tree shadows more difficult to distinguish from snowdrift, rock or shrub shadows (Supplementary Figs. 5 and 6). We suspect that many trees in the height class of 2–3 m were missed. These line segments, saved as .kml files, were imported into R (v.4.1.1)57 using the sf package58, where the length of each line segment was calculated and the coordinates of the shadow’s base were identified. The line segment lengths were used to estimate tree heights, and the coordinates were used in nearest-neighbour calculations and extractions of gridded data values. We estimated snow depth at 2.5–3 m because geolocated trees measured as ≤2.5 m in the field (see below) did not appear on imagery. We observed some trees taller than 2.5 m with no visible shadows on imagery, possibly buried in deeper snow or growing in shadows cast by terrain at the time of image capture. Thus, our estimates of adult populations may be underestimates, although there were also errors of commission where shrub shadows were mistakenly classified as spruce (see following).Digitizing and field validationTo estimate identification accuracy (Supplementary Information sections 1.2 and 1.3) among the 1,971 digitized shadows used for population reconstruction (enclosed by red rectangles in Supplementary Figs. 1–4), we visited 157 shadow locations first identified on imagery (8% of the 1,971) and located in the field with the built-in GNSS of late-model Apple iPhones (models 12 Pro Max, 12 Pro and second-generation SE) with positional accuracy in the open landscapes estimated at 3 m. At these 157 locations, 11 shadows were cast by very tall willows (7%). Of the 146 shadows confirmed as trees, 2 were dead (1%) and 1 had a recently broken top with green foliage on the ground. We added the length of the broken top to the standing height measured with a laser range-finder. Trees that were collinear in the solar azimuth at image capture contributed to errors of omission. The tree standing to solar azimuth obscured others as overlapping shadows fell in line, generating both errors of omission and an overestimate of the height of the first tree in the series. Six trees shadowed in three instances by what we identified on imagery as single shadows fell in this category. An additional three trees were missed during digitizing, also going unnoticed during QAQC, and were discovered in the field when matching shadows with trees. Supplementary Information section 1.3 provides details and a confusion matrix.In summary, 157 trees were expected from digitized shadows and 155 were found in the field. Applying the accuracy of the count overall suggests that 1,945 trees would better estimate the reconstructed population. Across the AOI, the total adult count of 5,988 shadows may represent 5,910 trees. Moreover, in so far as our estimates of ages based on tree heights are predictive, perhaps 2% of the ‘trees’ in our reconstruction are not a single tree casting a long shadow, but 2–3 younger, collinear trees. Thus, our estimate of past populations may be slightly biased to older trees, implying that the population growth rate may be slightly higher than estimated. However, the slightly fewer trees than shadows would suggest that the growth rate is lower. The relative size of these errors appears minor, and we did not incorporate them into the analysis, which seems to us robust and perhaps conservative in adult abundance estimates owing to image degradation with GEP super-overlays and other errors of omission. This study would have benefited from less image degradation using dedicated geographic information system (GIS) or image software. However, the low cost, simplicity and convenience of GEP was appealing for the large-scale digitizing.Returning from the field with individual tree data, R.J.D. displayed digitized shadow points together with field points on GEP, visually matching each field point to the nearest shadow, conditional on relative congruence between shadow size and tree height. This required care in clumps of trees with varying heights (example in Supplementary Information sections 1.2–1.3). The relative patterning of field points compared with shadows and the lengths of shadows compared with tree heights in these cases provided some measure of confidence in attribution.We made field expeditions to six study areas within the extent of the WV imagery we used for digitizing, three within the ‘simulated population area’ rectangle in Extended Data Fig. 1a (red rectangle in Supplementary Figs. 1–4) and three study areas further east (Extended Data Fig. 1c and Supplementary Fig. 2). Among-area variability was apparent in snow depth, terrain slope relative to the solar azimuth at the time of image capture and the solar-elevation angle itself because of the timing of image capture. The variability was identified, calculated and applied on the basis of geographic variability in the heights of trees casting shadows and from the slope and intercept of a mixed-model linear regression of field-measured height on digitized shadow length (see below).Field surveysWe validated species and heights of spruce casting shadows within the AOI along 403 km of ground transects. Our sampling did not appear spatially biased when compared with imagery as measured by proximity to a remote fixed-wing-aircraft landing site. Four field campaigns focused on three objectives in watersheds that were within or adjacent to the Noatak basin but did not have established treelines visible on WV growing season scenes: (1) to locate and document colonists at the geographic range boundary of white spruce; (2) to verify the locations of a sample of trees suggested by imagery in the AOI; and (3) to collect ecological measurements germane to white spruce range expansion. For adults (trees ≥2.5 m), datasets included height above ground (n = 340), diameter at breast height (DBH (~1.4 m); n = 296), CAG (n = 17), foliar nutrient content (n = 17), basal increment cores taken ≤20 cm above the ground (n = 140), tall shrub abundance within 5 m of sampled adults (n = 246), counts of juveniles within 5 m of sampled adults (n = 250), abundance class of cones (n = 339) and status of adults (live, n = 340; dead, n = 8). Of the dead adults, seven of eight were standing and largely without bark, with a median height of 4.1 m. The fallen dead tree was 6.2 m long with a DBH of 13.4 cm; all bark and limbs to fine branches remained. Only one dead adult, 4.1 m tall with a DBH of 4 cm, showed signs of decomposition with shelf fungus on the stem and decomposed limbs on the ground. Five juveniles ≥1.5 m tall had been stripped of their bark and all but their uppermost branches by apparently either porcupine (Erethizon dorsatum) or snowshoe hare (Lepus americanus). Anecdotally, we recorded other signs and possible causes of damage such as wind, bear (Ursus arctos), caribou (Rangifer tarandus) or struggling growth such as layering, stunted krummholz or clonal reproduction, although these growth forms were nearly totally absent.Field measurements for n = 770 juveniles located in the AOI and presented here included overall height, height above ground of bud scars representing 2015–2020 height (n = 302), damage and status. We used these measures to estimate age to increment core of adults (Supplementary Information section 2) and the RGR of juveniles (Supplementary Information section 3).Range expansion analysesDigitized established treelines (DETs) used here were downloaded as CTM_Treeline.kml from https://arcticdata.io/catalog/view/doi:10.18739/A2280506H. Ref. 34 describes drawing DETs on very high-resolution satellite imagery such as WV and Quick Bird. We clipped DETs to the four USGS HUC 10 watersheds within the HUC 8 Middle Kobuk subbasin and adjacent to the AOI (see ‘Environmental conditions’ below). The coordinates of the vertices for the clipped DETs provided the 3,366 locations of established treelines.We used the rdist.earth() function in the R package fields59 to identify the nearest neighbouring mapped adult and juvenile colonists in the AOI and DET vertices in adjacent Kobuk watersheds (Supplementary Information sections 1.8 and 1.9). Using the coordinates of nearest neighbours, we calculated differences in latitude as latitudinal displacement. Displacement north equalled the product of latitudinal displacement and 111.32 km, the distance between 67º and 68º N along 157.6891º W, which splits the AOI. Displacement in elevation was found by extracting from Interferometric Synthetic Aperture Radar (IFSAR) Alaska 5-m digital elevation models (DEMs) the elevation of DET vertices, mapped adults and mapped juveniles using the extract() function in the raster R package60 and then subtracting the elevation of the nearest neighbours from focal adults and juveniles. When geolocated adults or juveniles had estimated establishment years (see ‘Individual growth’ below), we calculated movement rates as the difference between the establishment year of an aged tree and the establishment year of the oldest tree sampled (1901, year of founding) as the denominator and displacement (difference in metres above sea level, kilometres or degrees of latitude) as the numerator (Supplementary Information sections 1.19–1.21). To time the progression of spruce away from DETs, we also binned establishment year by decade as decadal class, identifying within each decadal class the maximum displacement in kilometres north of and elevation in metres above (or below) nearest neighbours.Population growthFrom the 5,986 spruce shadow lengths within the AOI (Extended Data Fig. 1b and Supplementary Fig. 1) that we digitized from snow-covered scenes of DigitalGlobe WV imagery (Extended Data Table 1), we identified a sample of shadows stratified by length and cast by spruce that we located with GNSS-equipped late-model iPhones. We measured the height of n = 260 trees using a laser range-finder (LTI TruPulse 200) and/or a smartphone app (Arboreal Tree on iPhone 12 Pro and Pro Max with laser scanners) and collected n = 122 basal cores from individuals ≥2.5 m in height, then matched to shadows on imagery as described above (see ‘Digitizing and field validation’). Using the relationship between height and shadow length and the probability distribution of establishment year for the 122 cored trees identified within five height classes (Extended Data Fig. 2b), we simulated population growth within two contiguous sub-watersheds (the 135 km2 ‘simulated population area’in Extended Data Fig. 1a; western portion in Extended Data Fig. 2a; red rectangles in Supplementary Figs. 1–4; details in Supplementary Information section 4). These sub-watersheds contained n = 1,971 shadows cast on 26 March 2018. We treated these shadows as single spruce but recognize that they include as many as 138 willows (7%) and calculate an additional 118 (6%) spruce missed either by digitizing omission or by collinearity (Supplementary Information sections 1.2 and 1.3). Incorporating these errors together would not change the outcome of the simulations enough to change the doubling time of the population by more than a few percent.Estimates of tree height from shadow lengthOn a flat landscape covered uniformly in snow, the total height H of a tree equals snow depth S added to the product of shadow length L on the snow surface and the tangent of solar-elevation angle 𝛼, as H = S + Ltan(𝛼). However, because both the relative solar elevation and snow depth vary with terrain, we used a linear mixed-effects model (lmer() in the lme4 R package61) of height on shadow length (random factor of sample area with six levels), interpreting the fixed-effects intercept as the average snow depth (mean ± s.e. = 2.84 ± 0.14 m, t = 20.29) and the regression coefficient as the average tangent of solar elevation relative to the terrain slope (0.27 ± 0.04 m m−1, t = 6.96; details in Supplementary Information sections 4.1 and 4.2).Using these fixed-effects estimates and the random-effects covariance matrix, we applied Monte Carlo sampling to estimate the 1,971 heights with each run of the simulation, thereby propagating the error in height estimates. These 1,971 heights were then binned into five height classes with 0.5-m intervals from 4–5.5 m and with ≥1-m intervals from 3–4 m and 5.5–7 m (details in Supplementary Information sections 4.3 and 4.4). Height classes deduced from the shadow measurements were in some cases only 0.5 m in width. Because the mean snow depth (the intercept in the mixed-effects model) differed by more than this from one part of the study area to another (BobWoods, GaiaHill and BuffaloDrifts in Supplementary Information sections 4.1 and 4.2), this approach may have introduced systematic misclassification between locations. While applying a Monte Carlo model with coefficients drawn randomly using the mvrnorm() function from the MASS package in R with the random-effects covariance matrix was meant to alleviate this, we also ran the simulation with three uniform height classes with a wider interval (1.3-m width, for classes of 3–4.3 m, 4.3–5.6 m and 5.6–7 m).Estimating population-scale establishment yearWe estimated establishment years for each of the 1,971 trees (Supplementary Information sections 4.3 and 4.4). We did so by using the establishment yeardistributions by height class as Gaussian kernel densities for the 122 aged adults binned into the five height classes defined above (Extended Data Fig. 2b). Kernel density estimates were constructed using the function density() in R with options bw = “SJ” as the smoothing bandwidth, n = 107 as the number of consecutive establishment years, from = 1897 as the earliest year and to = 2004 as the latest year. For each of the 1,971 estimated heights binned into height classes, an establishment year was drawn (with replacement) from the corresponding kernel density distribution. We interpreted the total number of individuals in each establishment year as ‘recruitment by year’ into the population of survivors that we had digitized on the 2018 imagery. Sorting and cumulatively summing recruitment by year gave what we interpreted as population size (N) for each year (t) for trees that survived to 2018. Resampling in this manner for 1,000 runs, each time fitting exponential growth equation N(t) = N0ek(t – 1900) using nls() in R and then averaging the population RGR, provided population doubling time as ln(2) divided by mean k. The simulation was run again using three height classes, each of 1.3 m in width. The resulting mean doubling time was unchanged, but variability increased (Supplementary Information section 4.6).Individual growthCurrent annual growth and foliar chemistryIn autumn 2019, we collected current-year lateral branch tips on the west and east sides of each sampled spruce (n1 = 17 adult colonists and n2 = 457 adults at established treelines) at 1.4 m above the ground. Current annual branch growth was measured on 2–6 branches per spruce from the previous year’s bud scar to the tip of the branch. The number of samples varied, ensuring sufficient mass for foliar chemical analysis. Established treelines were sampled for adult foliage in 12 watersheds of the Noatak, Kobuk and Koyukuk river basins where we have ongoing experiments. At these sites, we used a replicated nested plot-based design (Extended Data Table 3). Colonist foliage sample locations (n = 8) in the upper Noatak basin were widespread across three watersheds. At each location, except the upper Noatak where 1–3 spruce per location were sampled, we sampled n = 5 white spruce separated by ≥10 m at a DBH of 8–12 cm. Needles from each branch tip were pooled by individual, dried for 48 h at 60 °C and weighed. Needles of individuals were pooled by treeline location after grinding to powder using a steel ball mill grinder (Mini-Beadbeater, Biospec Products) and subsampled for chemical analysis. Foliar N and 15N isotope were analysed for one subsample run on an Elemental Combustion Analyzer (Costech, 4010) coupled to an isotope ratio mass spectrometer (Delta Plus XP, Thermo Fisher Scientific) at the University of Alaska Anchorage Environment and Natural Resources Institute Stable Isotope Laboratory. Foliar P was measured for another subsample by the Pennsylvania State College Analytical Services Lab using the acid digestion method and analysed by inductively coupled plasma emission spectroscopy62.Juvenile RGRSeveral results presented here depend on juvenile vertical height growth during 2015–2020, which we assumed followed h(t) = h2015e(RGR t), where h(t) is height above ground for year t after 2015, h2015 is the height above ground in 2015 and RGR is the relative growth rate (Supplementary Information section 3). We used juvenile RGR in three contexts: (1) as a means of estimating establishment year in juveniles (Supplementary Information section 3.3); (2) as a metric of growth for comparison between colonist and established treeline juveniles (Supplementary Information section 6); and (3) to estimate the establishment year of cored trees (see second paragraph in ‘Dendrochronology’ below and Supplementary Information section 2).To estimate the RGR for each of 505 juveniles (n1 = 300 juveniles from m1 = 4 colonist populations and n2 = 205 juveniles from m2 = 14 established treelines; Extended Data Table 2), we measured the heights above ground (h) of the six uppermost bud scars in 2020, representing height increments in 2016–2020, the five consecutive years with the warmest mean daily July air temperature on record for Kotzebue. RGR in each juvenile was calculated as the regression slope of ln(h(t)) against t (mean R2 = 0.99 for 300 colonist regressions and 0.98 for 271 established treeline regressions; Supplementary Information section 3.4).To estimate the establishment year of juveniles, we used RGR to back-calculate T, the years required for an individual colonist to grow from 2 cm to h2015, as T = ln(h2015/2)/RGR. By subtracting T from 2020, we estimated the establishment year of each juvenile (Supplementary Information section 3.3).RGR values for colonist and established treeline juveniles (Extended Data Table 2) were compared using a linear mixed-effects model with field site (m = 24) as a random intercept, ln(RGR) as the dependent variable, ln(h2015) as a covariate to capture allometric growth and population (colonist or established treeline) as the fixed factor of interest (Supplementary Information section 6). Using the lmer() function of the lme4 package61 in R with REML = F, we found that the Akaike information criterion (AIC) for the interaction model was lower than that for the corresponding additive one (∆AIC = 22, likelihood ratio test χ2 = 24, degrees of freedom = 1, P 1 km beyond the established treeline, we recorded the location, age classes and presence of cones when possible. In watersheds of the uppermost Noatak basin and the Wulik basin, we also recorded both the total height of juveniles and the height above ground of the sixth bud scar from the tip to estimate RGR and so estimate age. We encountered three watersheds with tree island krummholz >1 km beyond the treeline but do not include these as colonist populations because clonal growth can be very old9,10,11,12,13,14,15,16,17,18,19. Of the 34 watersheds in which we encountered colonist populations >1 km beyond established treelines, 4 watersheds were located between 141° and 149.7° W (eastern Brooks Range), 21 watersheds were located between 149.7° and 156.3° W (central Brooks Range) and 9 watersheds were located between 156.3° and 163.3° W (western Brooks Range). Watersheds west of 150.5° W with colonists are shown in Fig. 1a.In 2021, R.J.D. led a field expedition to a small watershed in the Koyukuk basin (Arrigetch Creek, 67.439° N, 154.090° W). The watershed had been purposefully surveyed for juvenile white spruce above and beyond the treeline during 1978–1980 when seven juveniles 11–112 cm tall (six seedlings More
