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 area
The 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 expansion
Digitizing spruce shadows
We 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 validation
To 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 surveys
We 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 analyses
Digitized 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 growth
From 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 length
On 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 year
We 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 growth
Current annual growth and foliar chemistry
In 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 RGR
Several 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 < 0.0001). Applying lmer() with REML = T, we report the t score of the interaction between population and ln(h2015) as a test of the null hypothesis that mean RGR in colonist and established treeline juveniles was equal (Supplementary Information section 6.2).
Dendrochronology
We collected one basal increment core per tree for n = 140 trees at core height C, where 2 cm ≤ C ≤ 20 cm above ground (median, 5 cm). We selected trees with generally symmetric crowns in open locations with bark encircling the entirety of their mostly circular bole at core height. Increment cores were used to age individuals and to compare radial stem growth to mean daily July air temperature from Kotzebue during the period from 1990 to 2020. Cores were mounted, sanded and scanned at 1,200 d.p.i. We imported scans into CooRecorder (https://www.cybis.se/forfun/dendro/index.htm) for ring width measurements. We visually cross-dated cores63 and then checked cross-dating using COFECHA64 and CDendro (https://www.cybis.se/forfun/dendro/index.htm). In cores lacking pith, we used CooRecorder’s pith locator tool, which estimates distance and years to pith using the width and curvature of the innermost tree ring widths. Tree establishment year (Y) was estimated by subtracting the establishment age at core height from either the pith year or the year of the innermost tree ring minus the estimated years to pith. Growth age to core height was estimated using the height and age relationship from n = 83 colonists with 4 cm ≤ h2015 ≤ 20 cm and RGR calculated as described (see ‘Juvenile RGR’ above). For each of these 83 juveniles, we calculated T = ln(h2015/0.9)/RGR as years to 0.9 cm. The value T formed the dependent variable in a linear mixed model (m = 4 field sites as random factor) with ln(h2015) as the predictor variable. The fixed effects gave T = −7.31 + 10.46lnC, where C is core height (Supplementary Information section 2).
Raw tree-ring measurements were processed (Supplementary Information sections 7.3–7.5) to yield time series as BAI and as residuals from AR models65,66. BAI compensates for a possible decrease in ring growth increment with tree age by combining the size of the tree (radius, R) and its growth increment (∆r) as BAI = 2πR∆r. Although BAI is appropriate for direct analysis of growth, analyses of climate–growth relationships should also account for autocorrelation found in tree-ring series65,67. AR (‘pre-whitening’) results in series that approximate white noise with means of 1. The AR order of best fit was determined for each series by AIC score (AR order of 1 was most common). This method removes all but the high-frequency (interannual) variation in the series, which can then be compared to climate series. We also applied the AR approach to the July mean temperature data, but two common methods of selecting AR models (the Yule–Walker and maximum-likelihood methods) indicated that these data already approximated white noise (AR order of 0; Supplementary Information section 7.2). Thus, we retained the raw temperature data in subsequent analyses. All tree-ring detrending and standardizations were performed using the R package dplR68.
We assessed relationships between ring indices and July temperature using Pearson’s product-moment correlation of both ln(BAI) and AR with Kotzebue mean daily July air temperature from 1989 to 2019 using the cor.test() function in R with a two-tailed significance test (Supplementary Information section 7.5). We grouped the 140 increment cores as ‘juveniles’ <30 years old (n = 15) or ‘adults’ ≥30 years old (n = 125) for correlation with Kotzebue July temperature from 1989 to 2019. Adults must allocate structural carbohydrates to both growth and reproduction, whereas juveniles do not allocate to reproduction, suggesting that the strength and direction of radial growth responses to temperature may differ between age classes.
Environmental conditions
We used USGS HUC 10 and HUC 8 hydrological basins and watersheds to delineate our study AOI (https://water.usgs.gov/wsc/a_api/wbd/subbasin19/19050401.html; accessed from https://apps.nationalmap.gov/downloader/#/). The AOI was located in HUC 10 watersheds Kaluich, Cutler, Amakomanak and Imelyak of the HUC 8 Upper Noatak Subbasin. The DETs nearest the AOI were located in the Redstone, Miluet, Akillik and Hunt watersheds (HUC 10) of the Middle Kobuk Subbasin (HUC 8).
July temperature (Extended Data Figs. 3 and 6c) and November–March precipitation (Extended Data Fig. 6a) records from Kotzebue Airport (Fig. 1a) were accessed from the National Oceanic and Atmospheric Administration (https://www.ncdc.noaa.gov/cdo-web/search). July air temperatures (Extended Data Fig. 3) and October–April wind directions for the Kaluich (13 km west of the AOI; 758 m above sea level) and Imelyak (44 km east; 1,089 m above sea level) RAWS stations were accessed from the Desert Research Institute (https://raws.dri.edu/akF.html). Thirty-year (1980–2010) gridded (0.00833° resolution) mean July air temperature and November–March precipitation were accessed from the PRISM69 climate group (https://prism.oregonstate.edu/projects/alaska.php). We summed the PRISM November–March precipitation data (the water equivalent of snow for those months) by pixel. All correlations using these data sources applied Pearson’s product-moment correlation (r) with the cor.test() function in R and a two-tailed significance test. Strong, stable environmental lapse rates in July occur from Kotzebue to 1,200 m above sea level in the western Brooks Range (Extended Data Fig. 3 and Supplementary Information 7.2), making the instrumental temperature record there relevant to the AOI.
We obtained sea-ice cover in the Chukchi Sea from the National Snow and Ice Data Center (https://nsidc.org/arcticseaicenews/sea-ice-tools/). Time series of Chukchi Sea open water were derived from the file labelled ‘Monthly sea ice extent, by region’ (N_Sea_Ice_Index_Regional_Monthly_Data_G02135_v3.0.xlsx). We identified the overall maximum ice extent across years and months for the Chukchi Sea (maximum of 8.232 × 105 km2 in April 2006) and then subtracted October ice cover values from the maximum and defined the difference as ‘October Chukchi Sea open water’ (Supplementary Information section 8.3). We chose October as snow from that month is the first of the season that generally remains throughout the winter.
Bivariate climate envelopes in the AOI were based on the PRISM69 30-year mean July temperature and summed 30-year mean November–March precipitation extracted within the four HUC 10 watersheds in the HUC 8 Middle Kobuk Subbasin and the four HUC 10 watersheds of the Noatak AOI (Supplementary Information section 8.1). To construct a bivariate climate envelope for DETs, we extracted the Middle Kobuk Subbasin gridded PRISM data using the coordinates of the DET vertices. For the colonist population, we extracted the gridded PRISM data within the AOI using the coordinates of all mapped adult and juvenile colonists.
DEMs for estimating the elevation of geolocated adults, juveniles and DETs consisted of 5-m-resolution digital surface models collected using IFSAR Alaska and downloaded using USGS EarthExplorer (https://earthexplorer.usgs.gov/).
Regional extent of colonization
Range expansion requires a sequence of successful life history stages: seed dispersal, seedling establishment, sapling recruitment and adult sexual reproduction. Over 22° of longitude (141–163° W) of the Brooks Range, R.J.D. led field crews in search of ongoing range expansion by juveniles <1 m tall representing successful dispersal (‘seedlings’), juveniles >1 m tall representing successful establishment (‘saplings’) and individuals >2.5 m tall representing potential reproduction (‘adults’). The unit of sample was the watershed. Where one or more white spruce individuals (‘colonist populations’) were encountered >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 <1 m, one sapling ≥1 m) were located and mapped48. Our resurvey of upper Arrigetch Creek found 70 juveniles (52 seedlings, 18 saplings) and 19 adults. Near the mapped location of the two tallest juveniles in ref. 48, R.J.D. found a cone-bearing adult, as well as an additional eight cone-bearing adults elsewhere in the watershed, for a total of nine trees with cones among 19 adults up to 8 m tall. Four decades earlier, the tallest tree reported there had been 1.12 m tall.
Reporting summary
Further information on research design is available in the Nature Research Reporting Summary linked to this article.
Source: Ecology - nature.com
