Hindcasts of passerine density in arctic and subarctic Alaska suggest noncomplementary responses to shrub expansion by tundra- and shrub-adapted species

ABSTRACT

Shrub expansion is among the most conspicuous of the warming-related phenomena occurring within tundra regions, but its effects on vertebrates are not well understood. Among passerines, the habitats of shrub- and tundra-adapted species are expected to increase and decrease, respectively, by roughly equivalent amounts. However, shrub expansion is highly discontinuous, and these predictions are typically based on coarse-scale models of vegetation that permit limited spatial heterogeneity. We used a hindcasting approach to investigate spatial and interspecific heterogeneity in passerine densities through time, as a function of changes in Normalized Difference Vegetation Index (NDVI). We conducted line transect surveys of passerines between 2015 and 2019 at four sites in arctic and subarctic Alaska. Collecting the encounter locations permitted fine-scale assessment of covariates of density. We fit spatial distance sampling models to these data and then used the predicted relationship between passerine density and NDVI to hindcast changes in density since 1999. Our hindcasts suggested limited areal changes in the habitat of tundra-adapted species over the last twenty years while the habitat corresponding to tall shrub–adapted species has likely expanded. Thus, our results contrast with previous predictions due partially to the finer spatial scale of our inference, which allowed for portions of the landscape to resist shrub expansion.

KEYWORDS:

• Birds
• climate change
• density surface models
• NDVI
• spatial distance sampling

Introduction

Warming in the Subarctic and Arctic is amplified, inducing a cascade of impacts within tundra regions (Sturm et al. 2005; Post et al. 2009; Clegg and Hu 2010). Shrub expansion is among the most conspicuous of these phenomena and is evident across the circumpolar Arctic (Myers–Smith et al. 2011; Frost and Epstein 2014). Colonization by tall shrubs along riparian corridors and on hillslopes is the most visible example (Sturm, Racine, and Tape 2001; Tape, Sturm, and Racine 2006), but a general pattern of greening across tundra landscapes is also apparent (e.g., Verbyla 2008; Jia, Epstein, and Walker 2009; Pouliot, Latifovic, and Olthof 2009), driven largely by increases in the cover of both dwarf and tall shrubs (Fraser et al. 2014).

These changes in vegetation structure and composition are expected to affect the distribution of vertebrates (reviewed in Wheeler, Høye, and Svenning 2018). For example, the recent range expansions of moose (Alces alces gigas) and snowshoe hare (Lepus americanus) in Arctic Alaska appear to be linked to riparian shrub expansion in the region (Tape et al. 2016a, 2016b). In Subarctic Alaska, shrub-adapted passerines have exhibited rapid upslope shifts in distribution apparently tracking the expansion of shrubs into higher elevations (Mizel et al. 2016; Mizel, Schmidt, and McIntyre 2021). Such shifts in distribution are expected among shrub-adapted passerines, but the response of tundra-adapted species has received relatively little study. In one such study, Boelman et al. (2015) predicted a 20 to 60 percent increase in the shrub-dominated habitat of white-crowned sparrow (Zonotrichia leucophrys gambelii) and an equivalent decrease in the graminoid-dominated tundra habitat of Lapland longspur (Calcarius lapponicus) in the Brooks Range foothills of Arctic Alaska by 2050. However, these predictions were based on modeling vegetation change at a coarse spatial scale (a 4.5-km pixel) that fails to accommodate much of the spatial heterogeneity in the extent of shrub expansion (Boelman et al. 2015).

Though the signal of vegetation change in tundra regions is evident at coarse spatial scales (Goetz et al. 2005; Bunn and Goetz 2006; Jia, Epstein, and Walker 2006; Verbyla 2008), shrub expansion is highly discontinuous, occurring primarily along rock outcrops, hillslopes, floodplains, and areas of incised topography (Tape et al. 2012; Swanson 2013, 2015; Brodie et al. 2019). Soil characteristics (e.g., acidity, thaw depth, and drainage; Tape et al. 2012; Swanson 2015), herbivory (Olofsson et al. 2009; Christie et al. 2015), climate (Swanson 2015; Brodie et al. 2019), and disturbance history (Higuera et al. 2009) all mediate shrub distribution in tundra environments and, consequently, result in varied susceptibility of the landscape to shrub expansion. In Arctic Alaska, tussock tundra appears to be relatively resistant to shrub expansion and mesic sedge/Dryas tundra shows only modest evidence of recent increases in shrub density (Swanson 2013). Thus, failing to capture spatial heterogeneity in the extent of shrub expansion when forecasting passerine distribution may be biased toward predictions of complementary responses on the part of tundra- and shrub-adapted species.

Here, we hindcast passerine density as a function of changes in the annual maximum Normalized Difference Vegetation Index (NDVI; Tucker and Sellers 1986). We used data from ten species at four sites in Arctic and Subarctic Alaska, thereby permitting assessment of spatial and interspecific heterogeneity in changes in density. The passerine data were collected between 2015 and 2019 using distance sampling line transect methods and were unique in the sense that they comprise individual encounter locations rather than counts aggregated at the plot level, thereby allowing fine-scale assessment of spatial variation in density under a point process framework (Hedley and Buckland 2004; Mizel, Schmidt, and Lindberg 2018). We used these data together with the NDVI time series to model current distribution of passerine species and to make predictions about their density in 1999, approximately twenty years prior to our surveys. Annual maximum NDVI strongly correlates with shrub cover in tundra environments and is used widely for landscape-level assessment of changes in tundra vegetation structure (Jia, Epstein, and Walker 2003; Riedel, Epstein, and Walker 2005; Blok et al. 2011; Pattison et al. 2015).

Our hindcasts assumed that density in the past is a deterministic function of changes in NDVI. Therefore, our inference was to relative changes across species and sites due to changes in NDVI alone, which is analogous to predicting areal changes in potential habitat. That is, the predicted differences in bird density served as our unit of relative change in habitat. Given our focus on relative differences due to changes in breeding habitat alone, our objectives did not include predicting actual population trajectories that integrate numerous factors outside of breeding habitat structure. However, we expect that our hindcasts will help understand the extent to which current and future trends in the density of passerines may incorporate the effects of habitat change. We expected that our hindcasts would yield predictions of lower density for shrub-adapted species in 1999 compared to the present and we expected that species associated with tussock and alpine tundra would show limited predicted differences in density over this time period due to the relatively low susceptibility of these habitats to shrub expansion (Swanson 2013, 2015).

Methods

Study area

The study area is located within the U.S. National Park Service’s Arctic and Central Alaska Inventory and Monitoring Networks. We established four long-term monitoring sites targeting passerine assemblages breeding above treeline. Three sites are located at similar latitudes in the Arctic but are aligned along a longitudinal gradient extending from the central to western Brooks Range (Figure 1). The lone subarctic site is located more than 500 km to their south in the Alaska Range (Figure 1). Sites were located at access points including tundra airstrips in Noatak National Preserve (hereafter, arctic-west and arctic-central sites, respectively), the village of Anaktuvuk Pass within Gates of the Arctic National Park and Preserve (arctic-east), and along a gravel road in Denali National Park and Preserve (subarctic; Figure 1). We surveyed one site in each year from 2015 to 2019, resulting in only one of the four sites being resurveyed (arctic-east in 2015 and 2018). Because we were not specifically interested in interannual variation in density, we fit the same model to both years of data from the arctic-east site while allowing for site and year stochasticity through the incorporation of random effects (see below).

Figure 1. (a) Location of passerine monitoring sites in three national parks in Alaska including arctic-west (AW) and arctic-central (AC) in Noatak National Preserve, arctic-east (AE) in Gates of the Arctic National Park and Preserve, and subarctic (SU) in Denali National Park and Preserve. (b) An example of our route, plot, and grid configuration from the arctic-central site.

Mean elevations at the four study sites vary from 380 m at the arctic-west site to 1,140 m at the subarctic site (Supporting Information, Table S1). As a result of regional latitudinal and longitudinal temperature trends, the mean July temperatures at the four sites are similar and range from about 8.5°C to 11.5°C (PRISM Climate Group 2020; Supporting Information, Table S1). Based on ground sampling in the Arctic Inventory and Monitoring Network, shrub canopy volume (Canopy cover × Height of medium to tall shrubs) is typically low below 10°C mean July temperature and may exceed 1 m³/m² (volume per unit area) at 11°C or higher on suitable sites (Supporting Information, Figure S1).

Each site covers an elevational gradient from lowland tundra to alpine habitats including upland and riparian shrublands, herb- and low shrub–dominated wet tundra, mesic to dry Dryas and/or mixed dwarf shrub tundra, and barren slopes (Boggs et al. 2001; Jorgenson et al. 2009; Supporting Information, Tables S3 and S4). NDVI generally increases along a vegetation density gradient, from an average of 0.31 for barrens ecotypes to 0.54 for dwarf shrub (<0.2 m) ecotypes, 0.64 for herbaceous ecotypes, 0.68 for low shrub (0.2–1.5 m) ecotypes, and 0.77 for tall shrub (>1.5 m) ecotypes (Swanson and Neitlich 2016; after Viereck et al. 1992, Supporting Information, Table S5). NDVI varies widely within these basic structural classes due to the wide range in species and canopy covers that they encompass.

Though maximum NDVI is a coarse representation of the underlying vegetation gradient, it has fundamental relationships with ecological variables in tundra environments. This becomes apparent when we examine the shape of the NDVI–ecological variable relationship across studies. Specifically, we identified three studies in arctic tundra environments that fit exponential regressions with maximum NDVI as the independent variable and some ground-measured ecological variable as the response. The response variables were diverse and included arthropod biomass (Sweet et al. 2015), shrub aboveground biomass (Berner et al. 2018), and vegetation height (Bartsch et al. 2020). However, the shapes of these relationships are remarkably similar after standardizing the response to have mean 0 and standard deviation 1 (Figure 2). That is, the expected change in the standardized ecological response per unit of NDVI was similar across these studies. This implies that NDVI is capable of capturing variation in productivity and structure in tundra environments that is likely driven to a large extent by the density of deciduous shrubs.

Figure 2. Exponential relationships between maximum NDVI and ecological variables in arctic tundra. In the main plot, the relationships are shown for the same range of NDVI values (0.5–0.8) with the response (the ecological variable) standardized to have mean 0 and SD = 1. The insets show these relationships on their real scale within the range of NDVI values evaluated in each study and with line types corresponding to those shown in the main plot, although they overlap considerably. The three insets include (a) $binnedvegetationheight=0.38\times \mathrm{e}\mathrm{x}\mathrm{p}\left(8.1635\ast NDVI\right)$ (Bartsch et al. 2020), (b) $shrub\hspace{0.17em}aboveground\hspace{0.17em}biomass\hspace{0.17em}\left(AGB\right)=0.0012\times \mathrm{e}\mathrm{x}\mathrm{p}\left(8.89\ast NDVI\right)$ (Berner et al. 2018), and (c) $arthropod\hspace{0.17em}biomass=0.008\times \mathrm{e}\mathrm{x}\mathrm{p}\left(8.835\ast NDVI\right)$ (Sweet et al. 2015).

NDVI change

We derived NDVI rasters from 30-m resolution Landsat data (Landsats 5–8) in the U.S. Geological Survey’s Landsat Analysis-Ready Data (ARD) collection (Sauer and Dwyer 2018). The Landsat ARD collection contains atmospherically corrected surface reflectance data organized into a regular grid of 150 × 150 km tiles. We limited these data to all tiles with less than 40 percent cover of clouds/cloud shadow and those that were taken during the period of June through August, 1999–2019. Landsat coverage in our study area was sporadic until 1999, when multiple passes became available every summer.

We used the quality band within each Landsat ARD tile to designate as missing data all pixels classified with medium to high confidence as cloud, cloud shadow, snow, or ice. We corrected for slight differences in the band definitions of Landsat 8 using the method developed in Roy et al. (2016).

For each study area and year, we constructed a composite of the maximum NDVI values representing the peak growing season value. Peak NDVI correlates strongly with tundra vegetation structure, especially shrub density (Walker et al. 2005). Annual maximum values did not always represent the peak of the growing season due to the presence of cloud cover during the relatively short arctic summer and the scan-line corrector failure that affected Landsat 7 starting in 2003. Off-peak NDVI values reflect stages of leaf-out or senescence that do not properly represent the vegetation structure, which introduces error in peak NDVI values when ignored. That is, the raw, maximum NDVI values for each year contain a considerable amount of noise in relation to the peak greenness of deciduous shrubs, which is likely relatively static absent their incremental expansion. Because this was underlying process of interest in our hindcasts, we modeled peak NDVI using three-year composites to minimize this noise. Specifically, we aggregated the annual peak NDVI data from groups of three consecutive years (1999–2001, 2002–2005, …, 2017–2019), identifying each by the middle year, and computed the maximum of the three years to obtain a more reliable peak NDVI value. We then fit pixel-level linear regressions to the seven three-year composite values to obtain annual predictions of maximum NDVI using the “raster” package in R (Hijmans et al. 2015; R Core Team 2018). Our composite approach sacrificed some precision to reduce bias, which in turn resulted in more conservative estimates of hindcasted changes in bird density due to filtering changes in NDVI by significance level (see below).

Passerine surveys

We adopted a route-based, continuous effort survey design to maximize encounter rates of passerines. At the three roadless sites in the Arctic, we centered a grid comprising 600 × 600 m cells on the access point (Figure 1). We used this grid to establish a series of survey routes that radiated outwards from the access point (i.e., the basecamp/airstrip), permitting us to survey all accessible cells within 6.75 km (Figure 1). At the subarctic site, the grid was centered along the Denali park road and designated parking areas were used as access points for individual routes.

Routes were designed to be completed in a single day and were allowed to change direction in part to return observers close to their original starting point. In general, routes were located such that they intersected the centroid of each cell and then were deflected at an angle of 0° or 90° (i.e., followed a parallel or perpendicular direction; Figure 1). In some cases, we shifted routes off the central axis of the cell due to the presence of cliffs and high flow creeks. At each site, we surveyed twelve to seventeen routes an average of 3.6 times per season. Routes were fixed after their initial survey with one exception. During the pilot season (2015) at arctic-east, adjustments were made to the routes after their initial visit, necessitating survey-level indexing of the detection function in the analysis.

We conducted surveys of each route as if they were one continuous line transect, recording one observation per individual bird during a given route visit. Two-person survey teams included a data recorder who transcribed all encounter information spoken by the observer. Upon detecting a bird, we took a bearing and waypoint from our location along the line and used a laser range finder to estimate the distance to the bird. The encounter locations were projected in continuous space in ArcGIS using the waypoint, bearing, and distance data collected in the field. With one exception, the analysis was restricted to detections of singing males, which can be reasonably expected to represent the breeding population. For American pipit (Anthus rubescens), we pooled singing and calling detections based on similar detection functions estimated for both types of observations during exploratory modeling.

We delineated plots within each route by buffering the route using a species-specific truncation distance (based on removing ~5 percent of the furthest detections; Buckland et al. 2001) and clipping the buffered area at the intersection of the route and grid cell boundaries (see Figure 1). The purpose of separating the route-level survey strip into multiple plots was to accommodate variation in temperature over the course of a route visit (see below). We then overlaid a grid of 30-m pixels on each study site g= 1, …, G and assigned the observed encounter locations to the nearest pixel centroid leading to a data structure, which was a series of counts in pixel i= 1, …, I, in plot j = 1, …, J, and on visit t = 1, …, T (Mizel, Schmidt, and Lindberg 2018).

Analyses

Predicting the distributions of terrestrial vertebrate species in time and space is made difficult because vegetation trajectories are often relatively uncertain, themselves outcomes of multiple interacting processes. Specifically, the formation of new vegetation communities by shrub expansion at our sites is rate limited (by plant dispersal and growth rates) and limited at fine geographic scales by edaphic factors (Tape et al. 2012; Swanson 2015). Given the uncertainty associated with forward predictions of vegetation structure, we used a hindcasting approach to predict relative changes in passerine density across species and study areas due to the observed changes in maximum NDVI (hereafter, NDVI). At fine scales, NDVI explains much of the variation in passerine density in above-treeline areas (Mizel, Schmidt, and Lindberg 2018), in part due to the structural simplicity of the habitat and the distinctive reflectance patterns corresponding to barren ground, open tundra, and tall shrubs. On this basis, we modeled passerine density as a partial function of NDVI predicted for the survey year (2016–2019) from the pixel-level linear regressions. These models of current density then served as a baseline for hindcasting density with NDVI predicted for 1999 (the initial year in our NDVI timeseries).

Current passerine density

We fit spatial distance sampling models (Hedley and Buckland 2004) individually to each of ten species with sufficient data for assessing covariates of density at two or more of the sites. Our model yielded inference about density of the population that was present and available for detection. That is, we accommodated imperfect detection via the detection function, which describes the monotonic decrease in detectability with distance from the line. However, the repeated surveys were used to model variation in availability or song rate as nuisance parameters in a Poisson regression framework (Link and Sauer 1997; Oedekoven et al. 2013).

The focus was on estimating the intensity function for an inhomogenous Poisson point process of the individual encounter locations (Mizel, Schmidt, and Lindberg 2018):

(1) ${\lambda }_{ij}=exp\left({{\mathrm{Z}}_{ij}}^{\prime }{\beta }_g+{\epsilon }_k\right),$(1)

where ${\lambda }_{ij}$ is the expected number of individuals per pixel, the ${\beta }_g$ are the study site–specific vectors of $M\times 1$ coefficients, and ${{\mathrm{Z}}_{ij}}^{\prime }$ is the corresponding $\hspace{0.17em}M\times I\times J$ matrix of spatial covariates. The spatial covariates ${\beta }_g$ were made specific to each study site because elevations varied considerably across sites and to prevent one site from affecting the hindcasted predictions of the others.

The spatial covariates included the linear and quadratic effects of NDVI and elevation and the linear effect of solar radiation. The latter effect corresponded to the typical solar radiation (WH/m²) accumulated on 8 June (the average midpoint of our survey) given a generally clear sky. It was derived using the area solar radiation tool in ArcGIS and integrates variation in elevation, slope, and aspect, all variables that are expected to affect vegetation. Elevation was extracted from a 5-m digital terrain model (U.S. Geological Survey 2012). The NDVI values were those predicted for the survey year by pixel-level linear regressions of the NDVI timeseries.

We specified the normal random variables ${\epsilon }_k\sim N\left(0,{\sigma }_{\epsilon }^2\right)$ to accommodate repeated sampling within years and extra-Poisson variation. They corresponded to the 600-m grid cell or the 300 × 300 m subsection (as indexed by k) in which each pixel was located (see Figure 1). We used the coarser resolution grid (600 m) for sparsely distributed species where there was limited information about the distribution of these random variables (i.e., ${\sigma }_{\epsilon }$) at the finer scale (300 m). For the arctic-east site at which there was two years of data, we fit the same intercept and spatial effects but we specified a different ${\epsilon }_k$ for each grid cell–year combination, thereby accommodating interannual and spatial stochasticity.

For some study area–species combinations, we lacked sufficient detections for fitting a model with all five covariates. In such cases, we either dropped the solar radiation effect or fit an intercept-only model depending on sample sizes (see Supporting Information, Tables S6–S9). We did not estimate density for sites at which the species was not detected.

We described the observed encounter locations as arising from a spatial point process thinned through incomplete detection. The thinned point intensity ${\mu }_{ij}$ was the product of the expected point intensity and the detection function

(2) ${\mu }_{ij}={\lambda }_{ij}exp\left(-\frac{d_{ij}^2}{2{\sigma }^2}\right),$(2)

where d is the minimum distance between the route and each pixel and $\sigma$ is the scale parameter of a half-normal detection function (Mizel, Schmidt, and Lindberg 2018). We expected that the relatively dense vegetation at the arctic-west site would result in lower detectability for shrub-adapted bird species. Among this group of species and where sufficient data existed, we estimated a separate ${\sigma }_2$ for the arctic-west site while pooling the other three sites in estimating ${\sigma }_1$.

The observed pixel-level counts followed a multinomial distribution with cell probabilities ${\pi }_{ij}=\frac{{\mu }_{ij}}{\sum _i{\mu }_{ij}}$. The mean plot-level count was found by integrating the intensity function over the plot region using summation, ${\gamma }_j=\sum _i{\mu }_{ij}$ (Kéry and Royle 2016). We accommodated variation around the expected plot-level count in the model for the expected survey-level count ${\mathrm{\Lambda }}_{jt}$:

(3) ${\mathrm{\Lambda }}_{jt}=\mathrm{e}\mathrm{x}\mathrm{p}(log\left({\gamma }_j\right)+{{\mathrm{X}}_{jt}}^{\prime }\alpha ),$(3)

where $\alpha$ is a vector of $G\times 1$ coefficients and ${{\mathrm{X}}_{jt}}^{\prime }$ is the corresponding $\hspace{0.17em}G\times J\times T$ matrix of survey-level covariates. These covariates included the linear and quadratic effects of Julian date and the effect of temperature, assuming the presence of seasonal and diurnal variation in singing intensity (i.e., availability for detection; Mizel, Schmidt, and McIntyre 2021). The mean temperature for each plot visit was extracted from the weather station closest to each study site.

The observed survey-level counts were Poisson random variables, $n_{jt}\sim Poisson\left({\mathrm{\Lambda }}_{jt}\right)$, and the observed pixel-level counts were multinomial random variables conditional on the $n_{jt}$, $y_{ijt}|n_{jt},{\pi }_{ij}\sim Multinomial\left(n_{jt},{\pi }_{ij}\right)$.

Hindcasting passerine density

For each study site, we compared mean predicted density (males km⁻²) between 1999 and the survey year (2016–2019). To predict density in both years, we retained the full point intensity function (Equation (1)(1) ${\lambda }_{ij}=exp\left({{\mathrm{Z}}_{ij}}^{\prime }{\beta }_g+{\epsilon }_k\right),$(1) ) including all covariates and random effects. The only difference between the two predictions was that the one corresponding to 1999 was given the predicted NDVI values for 1999 in place of those predicted for the survey year with one exception. To minimize the potential for including any spurious changes in our hindcasts, we fixed NDVI at the survey year value for each pixel where we did not detect a significant trend from the linear regressions. We set the significance level at .10 for this filtering process.

Under this hindcasting framework, the effects of elevation and average solar radiation are assumed to be fixed over time. In predicting from Equation (1)(1) ${\lambda }_{ij}=exp\left({{\mathrm{Z}}_{ij}}^{\prime }{\beta }_g+{\epsilon }_k\right),$(1) alone, the estimates for both years (1999 and the survey year) reflect the survey-level covariates (survey temperature and date) held at their mean values. Because the baseline models for our hindcasts are estimated from one to two breeding seasons, their intercepts (i.e., the mean number of individuals per pixel) incorporate many factors affecting density, emphasizing our focus on relative changes in density due to NVDI alone.

Implementation

We fit all models in a Bayesian framework using Just Another Gibbs Sampler v4.3.0 (Plummer 2003) via the runjags package (Denwood and Plummer 2016) in R v3.4 (R Core Team 2018). We specified vague normal priors for all regression coefficients and half-Cauchy priors for the hierarchal standard deviation parameters (Gelman 2006). We estimated the posterior distribution from four Markov chain Monte Carlo chains run for 20,000 to 60,000 iterations. We assessed convergence visually and using the Gelman-Rubin diagnostic (Brooks and Gelman 1998). We scaled continuous covariates (mean = 0, SD = 1) to improve convergence properties. Inferences about changes in passerine density were based on whether 95 percent highest posterior density credible intervals overlapped 0.

Results

NDVI response

Across the four sites, changes in NDVI occurred over different parts of its range and in different directions, implying a varied set of vegetative transitions and marked spatial heterogeneity in the extent of shrub expansion. The subarctic and arctic-west sites exhibited increases in the proportions of their area with high NDVI (0.7–0.9), coupled with decreases in their proportions in the intermediate range (0.5–0.7; Figures 3a, 3d),implying increasing tall shrub cover. These transitions to higher NDVI classes were concentrated on mountain slopes, but increasing NDVI was also observed on dry sites such as stream gravel bars and convex bedrock hills that had low initial NDVI (0.1–0.5; Figures 3a, 3d).

Figure 3. Change in annual maximum NDVI between 1999 and the survey year in the (a) arctic-west, (b) arctic-central, (c), arctic-east, and (d) subarctic study areas.

The arctic-east site showed increases in the proportion of its area with intermediate/high NDVI (0.6–0.8), coupled with roughly equivalent decreases in the proportion of its area in the range of 0.3 to 0.6, implying increases in low shrub and herbaceous cover and modest tall shrub expansion (Figure 3c). Increases in high NDVI were largely restricted to riparian zones including gravel bars and water tracks.

In contrast, the arctic-central site showed decreases in the area of low (0.2–0.5) and intermediate/high NDVI (0.7–0.8) in concert with increases in intermediate NDVI (0.6–0.7), implying a net reduction in the areas of both unvegetated and shrub cover types accompanied by an increase in herbaceous cover (Figure 3b). Pixels that showed negative trends from initially high NDVI were concentrated in wet footslope areas where inspection of high-resolution aerial photographs reveals a loss of shrub cover, perhaps due to increasing wetness (Supporting Information, Figure S2). These losses in the intermediate/high NDVI range that occurred on wet sites were partially offset by gains on gravel bars and hillslopes. Increases in intermediate NDVI were observed at a wide range of mostly well-drained sites.

Spatial variation in passerine density

Covariates of distribution (i.e., NDVI, elevation, and solar radiation) and availability (i.e., date and temperature) were well supported across species and sites (Supporting Information, Tables S6–S10). Species were separated into broad habitat associations as inferred by NDVI (Figure 4). American pipit occupied alpine tundra habitats (low NDVI), Lapland longspur used graminoid-dominated tussock tundra (intermediate NDVI), and savannah sparrow (Passerculus sandwichensis) used environments characterized by low shrubs interspersed with herbs (intermediate/high NDVI). The remaining seven species used high NDVI environments and, to a lesser extent, low and/or intermediate NDVI environments, reflecting an association with open shrublands (Figure 4b). The shapes of these curves varied, indicating differences in thresholds for minimum shrub cover. The two most gradually sloping curves, belonging to white-crowned sparrow and American robin (Turdus migratorius), showed nonnegligible density in low NDVI environments and, therefore, some association with sparsely vegetated cover types. The curve for American tree sparrow (Spizella arborea) had an intermediate slope. The steepest curves, corresponding to golden-crowned sparrow (Zonotrichia atricapilla), Wilson’s warbler (Cardellina pusilla), and fox sparrow (Passerella iliaca), suggest specialized use of more extensive tall shrub thickets typically found in riparian zones (Figure 4b).

Figure 4. Density–NDVI relationships at the arctic-central site for (a) tundra and low shrub–adapted passerines and (b) tall shrub–adapted passerines. Curves represent the expected change in density per unit change of maximum NDVI, while holding the intercept and all other covariates at 0; that is, they are predictions of the exponentiated coefficients for linear and quadratic NDVI. The intercept has been removed to allow plotting the effects from all species on the same scale. Though the arctic-central site was used for display purposes, the shape of these relationships was similar across sites. Species codes in (a): American pipit (AMPI), Lapland longspur (LALO), and Savannah sparrow (SAVS). Species codes in (b): white-crowned sparrow (WCSP), American robin (AMRO), American tree sparrow (ATSP), golden-crowned sparrow (GCSP), Wilson’s warbler (WIWA), and fox sparrow (FOSP).

Passerine hindcasts

Densities of tundra-adapted species (American pipit and Lapland longspur) were predicted to be stable from 1999 to present with the exception of a small predicted increase in Lapland longspur density at the arctic-central site (Figure 5). Nearly all tall shrub–adapted species were predicted to have increased density at subarctic and arctic-west sites, and the size of these predicted changes tended to increase according to the strength of their association with dense tall shrubs or high NDVI (Figures 5 and6a, 6d); Supporting Information, Figures S3–S9). At the arctic-east site, predicted increases in the density of tall shrub–adapted species were less pronounced due to less change in high NDVI (Figure 3c) and/or imprecise estimation of change resulting from the low densities of these species at this site (Figure 5c; Supporting Information, Figures S10–S13). Nevertheless, three of five tall shrub–adapted species were predicted to have increased density at the arctic-east site (Figures 5 and 6c).

Figure 5. Relative change in density from 1999 to the survey year for passerine species at four sites in Alaska including (a) arctic-west (1999–2017), (b) arctic-central (1999–2016), (c) arctic-east (1999–2018), and (d) subarctic (1999–2019). Species are numbered according to the strength of their association with increasing NDVI; that is, they are ordered according to the average NDVI value (across sites) at which 50 percent of their total expected abundance was reached. Species included (1) American pipit, (2) Lapland longspur, (3) Savannah sparrow, (4) white-crowned sparrow, (5) American robin, (6), American tree sparrow, (7) golden-crowned sparrow, (8) orange-crowned warbler, (9) Wilson’s warbler, and (10) fox sparrow. Missing estimates for some species–study site combinations indicate that the species was absent from the site or had an intercept-only model. Posterior means and 95 percent highest posterior density credible intervals are shown in relation to a horizontal line at 1, which indicates no change.

Figure 6. Difference in point intensity (number of males/30-m pixel) for American tree sparrow between 1999 and the survey year within a 300-m sampling strip at four sites in Alaska: (a) arctic-west (1999–2017), (b) arctic-central (1999–2016), (c) arctic-east (1999–2018), and (d) subarctic (1999–2019).

The arctic-central site differed from the overall pattern in that all tall shrub–adapted species were predicted to have either stable or decreasing density, whereas the tundra- and low shrub–adapted species were predicted to have stable or increasing density (Figure 5b; Supporting Information, Figures S10–S16). This individualistic pattern likely reflected the thinning or loss of tall shrub thickets rather than their expansion.

Savannah sparrow, our lone species associated with low stature shrubs, was predicted to have stable mean density at two sites (Figures 5a, 5c) that masked fine-scale spatial variation in the trajectory of NDVI (Figure 7). Specifically, its association with intermediate/high NDVI (Figure 4a) coincided with the range in which considerable change occurred. At the arctic-west site, these changes appear to have corresponded to a general pattern of greening that resulted in offsetting areas transitioning in and out of this species’ optimal NDVI range (Figures 3a and 7a). At the arctic-central site, increasing NDVI on mostly well-drained sites and decreasing NDVI on poorly drained sites also resulted in predictions of offsetting areas of increasing and decreasing density (Figure 7b).

Figure 7. Difference in point intensity (number of males/30-m pixel) for Savannah sparrow between 1999 and the survey year within a 300-m sampling strip at two sites in Arctic Alaska: (a) arctic-west (1999–2017) and (b) arctic-central (1999–2016).

Discussion

By virtue of assessing both NDVI change and variation in passerine density at fine spatial scales, we permitted noncomplementary predictions of habitat change among shrub- and tundra-adapted passerines. Tall shrub–adapted species were predicted to have increased habitat area at arctic-west and subarctic sites, comparatively modest increases at the arctic-east site, and decreasing or stable habitat area at the arctic-central site. Conversely, the habitat area corresponding to tundra-adapted species was predicted to have changed little over the last twenty years. Our hindcasts contrast with predictions made by Boelman et al. (2015), who expected equivalent areal decreases and increases in the habitat of the tundra-adapted Lapland longspur and the tall shrub–adapted white-crowned sparrow, respectively. The disagreement in these predictions appears to be due in some part to the finer spatial scale of our inference, which allowed for portions of the landscape to resist shrub colonization presumably as a function of soil limitations (Tape et al. 2012; Swanson 2015). The tussock tundra habitat used by Lapland longspur has soil characteristics (e.g., relatively high moisture and acidity) that make them unfavorable for the formation of tall shrub canopies even under warmer conditions (Tape et al. 2012; Swanson 2015).

Our predictions are consistent with a growing body of evidence of individualistic responses to climate change. Although studies of range dynamics predominantly report a pattern of poleward or upslope movement in response to climate warming (Parmesan 2006; Hitch and Leberg 2007; Maggini et al. 2011), heterogeneity in the response across species and populations is also apparent, with numerous examples of unexpected responses including range stability, isotropic change (e.g., due to habitat fragmentation), and shifts that are opposite of the predicted direction (e.g., Lenoir et al. 2010; Tingley et al. 2012; Mizel et al. 2016). In addition to predicting relatively stable habitat area for the two tundra-adapted species that we considered, our lone species associated with low shrubs, savannah sparrow, showed limited evidence of areal changes in habitat at two of our sites. However, this was an outcome of offsetting areas transitioning in and out of this species’ optimal NDVI range over the study period. Thus, savannah sparrow is indicative of how coarse models of distribution may mask change at fine spatial scales, particularly in species with unimodal or humped responses to environmental variables.

In assuming that density in the past is a deterministic function of changes in NDVI alone, our approach was analogous to predicting areal changes in potential habitat. However, there are limitations to our inferences about habitat change due to the simplicity of our approach, including the fact that it assumes that the NDVI–density relationship is fixed over time. It also ignores potential changes on scales larger than a 30-m pixel that may have impacts on habitat suitability (Henden et al. 2013; Boelman et al. 2016). Specifically, shrub-adapted species showed use of open shrublands (i.e., intermediate/high NDVI), and continued infilling could result in suboptimal levels of canopy closure at patch scales or larger (Mizel et al. 2016), leading to local population declines in the reverse of our hindcasts. Along these lines, Mizel et al. (2016, 2021) documented upward shifts in the elevational distributions of nearly all common shrub-adapted species from 1995 to 2019 at the subarctic site using a timeseries of roadside counts, suggesting potential tracking of open shrub canopies over time. In addition, Thompson et al. (2016) found that some tall shrub–adapted species in western Alaska declined in abundance above a threshold level of average shrub height and/or shrub cover within 250-m radius plots. For some shrub-adapted species, openings in the shrub canopy provide herbaceous cover that aids in nest concealment and provides protection against severe weather (Boelman et al. 2015). Some shrub tundra species also have foraging tactics adapted for open shrub canopies (Batova 2011), and others show associations with structurally diverse shrub communities (Boelman et al. 2016).

Although tall shrub expansion was apparent to varying extents at all the study areas, other vegetative transitions were also implied by the observed changes in NDVI. The arctic-central site appears to be relatively resistant to shrub expansion and instead showed decreases in the cover of high NDVI that potentially corresponded to the thinning of shrubs in areas of wet tundra. Simultaneously, increases in intermediate NDVI occurred on sparsely vegetated dry sites, suggesting increasing cover of a mixture of tundra life forms (i.e., graminoids, forbs, dwarf shrubs, and erect shrubs). This is consistent with the observed pattern of increasing graminoid and dwarf shrub cover in the Arctic (Jia, Epstein, and Walker 2006; Joly, Cole, and Jandt 2007; Munger et al. 2008; Tape et al. 2012; Fraser et al. 2014). Therefore, the low susceptibility of tussock tundra habitats to shrub encroachment coupled with potential increases in low vegetation on barrens suggests that areal changes in the potential habitat of tundra-adapted passerines within our study area are unlikely to simply be the reverse of those corresponding to shrub-adapted species. Our study also emphasizes uncertainty about the future population trajectories and distribution of shrub-adapted species. The association of some of these species with open shrublands suggests vulnerability to continued infilling of shrub canopies and, potentially, the transitory nature of their density trajectories at some sites. In addition, our study stresses the importance of framing inference in terms of the individual encounter in space rather than aggregating counts at the plot level. This approach more fully leverages the spatial information underlying the individual encounters in explaining distribution (Mizel, Schmidt, and Lindberg 2018) but is also critical for accommodating spatial heterogeneity in the response of vegetation in predictive models.

Author contributions

JDM conceived and designed the experiments. JDM and DKS analyzed the data and wrote the article.

Acknowledgments

We thank the Nunamiut and Arctic Slope Regional Corporations for access to their lands. We thank the technicians who conducted surveys over the course of this project, particularly Jared Hughey, Sarah Swanson, and Shelli Swanson. We thank Carol McIntyre for help with logistics, Ken Hill and Pamela Sousanes for help acquiring weather data, and Josh Schmidt, Eric Wald, and two anonymous reviewers for comments on drafts of this article.

Disclosure statement

No potential conflict of interest was reported by the authors.

Supplementary material

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

Additional information

Funding

This study was funded by the U.S. Department of the Interior, National Park Service, Inventory and Monitoring Program in the Central Alaska and Arctic Inventory and Monitoring Networks.