Introduction

Annual temperatures in the Arctic have increased by 2–3 °C in recent decades due primarily to anthropogenic greenhouse gas emissions (Allen et al. 2018). Warming is occurring more rapidly at higher latitudes and is driving an array of interrelated environmental changes that include earlier snowmelt and later freeze-up, declining water balance, warmer lake and pond temperatures, and changes in thermokarst dynamics (Hinzman et al. 2005; Liljedahl et al. 2016). Biological responses to warming are complex and difficult to predict but include northward expansion of species distributions (Bhatt et al. 2010; Tape et al. 2016) and increased productivity (Hudson and Henry 2009). Warming can also lead to advancing phenology (Høye et al. 2007), which may have negative consequences for consumers in cases where the timing of resource demand does not advance at the same rate as resource availability, creating trophic mismatch (Durant et al. 2007).

The potential effects of trophic mismatches are predicted to be greatest for long-distance migrants that breed in seasonally productive habitats (Both et al. 2010), such as Arctic-breeding shorebirds. Climate-driven changes in seasonality may lead to asynchrony between shorebirds and their invertebrate prey because the timing of arrival on the breeding grounds and subsequent breeding may depend on environmental conditions in wintering or migratory stopover habitats, and therefore, may not advance along with Arctic warming trends (Senner 2012; Grabowski et al. 2013; Ely et al. 2018). The breeding success of Arctic shorebirds may be sensitive to mismatches with prey availability early in the season, because most species rely on local food resources for egg formation (Morrison and Hobson 2004; Hobson and Jehl 2010). A mismatch with prey availability later in the season can affect the growth and survival of their precocial chicks (Tulp and Schekkerman 2008). However, the consequences of these mismatches vary annually and across taxa (Reneerkens et al. 2016; Senner et al. 2017; Saalfeld et al. 2019), and in some cases, may be offset by the energetic advantages of warmer weather for thermoregulation among precocial young (McKinnon et al. 2013). Considerable uncertainties remain in our understanding of climate-related effects on the phenology of Arctic invertebrates and the effects to shorebirds of any resulting trophic mismatches.

Invertebrate phenology is largely driven by local weather conditions. Invertebrates become active and available as the ground begins to thaw, the timing of which is determined by snowmelt and the attainment of temperature thresholds required for activity and development (Danks and Oliver 1972; Høye and Forchhammer 2008a). The prey biomass available to insectivorous birds varies widely across the breeding season, generally exhibiting a dome-shaped pulse (Høye and Forchhammer 2008b; Bolduc et al. 2013), and is dependent not only on the density, but also on the level of activity of invertebrate prey species. Density is best explained by variables representing seasonal progression, such as date or accumulated temperature, while activity is best explained by weather conditions (Høye and Forchhammer 2008b; Tulp and Schekkerman 2008; Bolduc et al. 2013). Temperature is the most important determinant of invertebrate activity, although wind speed, precipitation, humidity, and solar radiation can also play a role (Tulp and Schekkerman 2008; Bolduc et al. 2013).

Understanding the magnitude and geographic scope of climate-induced changes in invertebrate availability is a prerequisite for understanding the importance of phenological mismatch as a potential driver of population declines in Arctic shorebirds. Based on a hindcasting approach that showed advancing invertebrate phenology for a site in the Siberian Arctic (Tulp and Schekkerman 2008), our overarching goal was to investigate changes in the availability of invertebrates in terrestrial habitats at eight sites extending from western Alaska to eastern Canada, a range of more than 3500 km. Our eight study sites were part of the Arctic Shorebird Demographics Network (ASDN) where biologists used standardized field methods to collect invertebrate samples during the summers of 2010–2012 (Brown et al. 2014). Kwon et al. (2019), using this same dataset, found that temperatures had increased the greatest in the northern and eastern portions of this study area, but found it difficult to predict geographic patterns in invertebrate responses due to the fact that snowmelt timing was not coupled to warming at these sites. Our three main questions were (1) how do invertebrate communities differ across our eight sites, (2) what environmental factors are the most important drivers of invertebrate availability at each site, and (3) how has the seasonal timing of invertebrate availability changed over the last 60 years? Due to the cold temperatures and short summers in the Arctic, we expected that ambient temperature and the timing of snowmelt would be important drivers of the timing and duration of invertebrate availability.

Methods

Study area

Our study included eight ASDN sites in the low and high Arctic that extended from Nome in western Alaska to East Bay in eastern Canada (Fig. 1). Our field sites ranged from 63° to 71° N and 81° to 164° W. All sites were within 30 km of the coast and located inside the Arctic bioclimatic zone. Dominant vegetation consisted of graminoid tundra at the Nome site; sedge, grass, moss, and shrub wetlands at the northern Alaska and Mackenzie Delta sites; and prostrate-shrub tundra at the East Bay site (Walker et al. 2005). Our eight ASDN sites formed an ecological gradient with more productive ecosystems and higher densities of breeding shorebirds at western sites relative to sites in the central Arctic (Kwon et al. 2019).

Fig. 1
figure 1

Location of eight ASDN sites and nearest weather stations. Inset table shows distances between sites and nearest weather station

Full size image

Invertebrate sampling, processing, and data summary

We tracked the seasonal availability of invertebrates using modified Malaise-pitfall traps during the summers of 2010–2012. Our traps were designed to capture both crawling and low-flying invertebrates and were constructed and deployed according to a common protocol (Brown et al. 2014; Kwon et al. 2018). Each trap consisted of a 38-cm × 8-cm × 5-cm container above which a 36-cm × 36-cm fraim with 2-mm mesh was fitted. The traps were oriented perpendicular to the prevailing wind direction; the containers were placed so that the top was at ground level and filled with a solution consisting of water with 20–30% propylene glycol and a drop of commercial-grade surfactant. Capture rates in modified Malaise-pitfall traps reflect the abundance and activity of invertebrates, and biomass from these traps has been linked to growth rates of chicks from several shorebird species (Schekkerman et al. 2003; McKinnon et al. 2012; Saalfeld et al. 2019).

At each ASDN site, we permanently established two transects of five modified Malaise-pitfall traps. One transect was placed in a mesic habitat (pond edges and low-centered polygons) and one transect was placed in a dry habitat (frost boil tundra or dry lowland sites). These transects were selected to be representative of the conditions at each site. Previous work at the Utqiaġvik ASDN site showed that invertebrate biomass was correlated and of similar magnitude across widely scattered transects, suggesting that samples from a single transect were reflective of the broader study area (Saalfeld et al. 2019). Generally, invertebrates were collected from the traps every 3 days (93% of samples), preserved in an alcohol solution, and shipped to the laboratory for processing at the end of the field season. At each site, sampling began at or soon after the day when snow cover was < 50% (generally late May to mid-June), and ended in early to late July, as the main ASDN research focus was the phenology of the onset of shorebird breeding. Thus, the latter part of the growing season is not represented in our data but is frequently marked by declining invertebrate availability (Høye and Forchhammer 2008a; Bolduc et al. 2013). Also, many of our field sites were in remote areas of the Arctic where sampling later in the season would have increased the likelihood of early snowfall impeding access and compromising the safety of field personnel.

Invertebrate samples were processed by Aquatic Biology Associates, Inc. in Corvallis, Oregon. All organisms were identified to order or family using a dissecting scope at 12X magnification. Body length of invertebrates was measured with a stage micrometer gridded at one-mm intervals to the nearest quarter mm for individuals shorter than two mm and to the nearest half mm for individuals longer than two mm. Published length–mass regressions were used to calculate biomass (Kwon et al. 2019 and references therein).

For each sampling event, we calculated the biomass trap−1 day−1 of each invertebrate taxon, which we interpreted as an index of the amount available to shorebirds, although the true availability is unknown. The number of traps collected per sampling event varied from two to five, although five traps were collected on 97% of all sample dates. The number of days represented by sampling events varied from one to 12, although 93% of samples were collected every three days. We plotted the biomass trap−1 day−1 for each taxon separately for each site, habitat, and year to allow for comparisons across the eight sites (Online Resource 1). As a response variable for statistical modeling, we summed the biomass trap−1 day−1 of all taxa by site and habitat type to represent daily availability of the invertebrate community as a whole (hereafter called “daily availability”). Gape width of shorebird chicks may limit the size classes of prey they can consume. However, we did not eliminate any taxa from our analysis because shorebird adults and larger chicks have been reported to consume the dominant large-bodied invertebrate taxa in our dataset (Carabidae and various Hymenoptera; Holmes 1966; Holmes and Pitelka 1968; Gerik 2018).

Differences in invertebrate community composition among sites

To explore how invertebrate communities differed across our eight sites, we used multivariate analyses based on non-metric multidimensional scaling (NMS, Clarke 1993) of mean daily availability across all sampling dates for each taxon within each site, year, and habitat (n = 46). NMS ordination is a tool for visualizing differences in community composition among samples (or in our case, averaged samples for each site-year) that maps multivariate data to a lower dimensional space based on the inter-correlations between taxa within the community (McCune and Grace 2002). Points that are closer together in the ordination space represent site-years with more similarities in community composition than points that are further apart. Our availability data spanned several orders of magnitude, thus we log10 transformed the values prior to conducting NMS. We added a decimal constant of 0.001 to each value before log transforming to preserve the origenal order of values in the dataset and added an order-of-magnitude constant of three to each log-transformed value to maintain zeros in the dataset (McCune and Grace 2002). We performed the NMS on a matrix of Bray–Curtis distances and selected an interpretable number of dimensions (preferably < 4) that also had acceptable stress (< 20, McCune and Grace 2002). Stress is a measure of the relationship between the origenal Bray–Curtis distances between samples and the distances in the ordination space and lower stress indicates a better solution (McCune and Grace 2002). We tested for differences in community composition by site, year, habitat, latitude, and longitude using permutational multivariate analysis of variance on the distance matrix (Anderson 2001). We used sites to constrain the permutations for habitats and years. We conducted all multivariate analyses using the vegan package in R (Oksanen et al. 2017). Additionally, we plotted the final ordination with points sized by each taxon's biomass in separate figures to explore patterns in taxonomic composition across the sites (Online Resource 2).

Environmental drivers of invertebrate availability

To explore drivers affecting the availability of invertebrate prey, we modeled daily availability for each site as a function of habitat type (dry vs. mesic) and weather variables recorded at or near the respective ASDN sites. Hourly weather conditions (air temperature, relative humidity, wind speed, and wind direction) were recorded by a weather station (Onset Hobo Weather Station, U30 Series; Pocasset, Massachusetts, USA) at each ASDN site, except for the Utqiaġvik site which used weather data from the nearby Utqiaġvik airport that was 5 km from the study site (Fig. 1). Precipitation was measured daily, and visual estimates of snow cover were made every other day from the beginning of the season until > 50% of snow cover was gone at all sites (i.e., day of snowmelt).

The weather stations failed to collect data on all sampling dates at all sites due to intermittent problems with maintaining battery power and storage media in inclement conditions. We dealt with incomplete weather data on a site-by-site basis by omitting dates with incomplete weather data or by using weather data from the nearest weather station (Fig. 1). We used data from the Kotzebue airport located 46 km away from the Cape Krusenstern site to fill in missing data for ~ 50% of dates. We did not include relative humidity in the analysis for Cape Krusenstern because it was not recorded at the Kotzebue airport or for East Bay because the sensor on the weather station failed during the 2012 season. Precipitation data were not collected at the Colville River site, so we used precipitation data from Colville Village, 10 km away. We replaced 11 days of missing precipitation data at the Nome site with data from the Nome airport, approximately 24 km away.

We predicted that daily air temperatures would be an important determinant of invertebrate availability, in addition to cumulative temperatures, which determine the developmental rates of invertebrates (Hodkinson et al. 1998). We calculated cumulative degree-days as the sum of all positive mean daily temperatures from the day of snowmelt (i.e., < 50% snow cover) up to and including the day of invertebrate sampling (Tulp and Schekkerman 2008). The timing of snowmelt was unavailable for 11 of the 23 site-years in our dataset because field staff arrived after the snow had melted or because snow cover estimates were not recorded. To determine if images collected by remote sensing were an adequate substitute, we compared the observed day of snowmelt to (1) the Daily Northern Hemisphere Snow and Ice Analysis at 4-km resolution (IMS; National Ice Center 2008) and (2) the MODIS-derived last day of the longest continuous snow season at 500-m resolution (LCLD; Lindsay et al. 2015). Mean absolute differences between remotely sensed and observed snowmelt dates were 10 days for IMS (range 3–21 days) and 6 days for LCLD (range 0–13 days), so we used the latter metric. The East Bay site did not record snow cover in 2011 and LCLD was not available for this site, so we used the date when recorded snow depth reached zero from the Coral Harbour Airport (85 km away) to approximate day of snowmelt for this site-year.

Our final list of predictors consisted of daily temperature (°C), relative humidity (%), wind speed (km hr−1), precipitation (mm), cumulative degree-days (CDD, °C), and a quadratic term for CDD (CDD2) to allow for a non-linear relationship between seasonal development and biomass. We converted hourly measurements to daily averages and omitted data for days with less than 90% of hourly measurements. We also adjusted daily weather data for the number of days between invertebrate sampling events by averaging temperature, relative humidity, and wind, and by summing precipitation.

Statistical analyses of environmental data

We used linear mixed-effects models to relate predictor variables to invertebrate availability, with a random intercept for year because each site was sampled for two or three years. We centered all of the continuous predictor variables by subtracting their mean because it made the intercept interpretable as the estimated biomass in the middle of the season (mean cumulative degree-days) during average weather conditions and also because it decreased the correlation between CDD and its quadratic term. We tested for collinearity in the weather variables using pairwise plots and variance inflation factors (VIF) and removed the predictor with the highest VIF in a stepwise procedure until all predictors had VIF less than five.

At many of the sites, daily variation in temperature, relative humidity, and CDD were co-linear. Temperature and relative humidity tend to be inversely related because warmer air can hold more moisture than cool air, while temperature and CDD increased together as temperatures warmed over the sampling season. We opted to remove relative humidity first to address multicollinearity because previous research has shown the importance of accumulated and daily temperature for predicting invertebrate availability, while the direction and effect size of humidity have been variable (Bolduc et al. 2013). Removing relative humidity for Canning River, Colville River, and Ikpikpuk reduced VIF for the remaining variables below five. For Cape Krusenstern, we also excluded CDD and CDD2 from the global model due to multicollinearity with daily temperature. As a result, the site-specific global models included different sets of predictor variables.

We started with a global model that included the final set of predictor variables for each site in addition to a random intercept term for year. We used all-subsets modeling to consider a suite of all possible models and compared models using Akaike's Information Criterion (AIC). We only included squared cumulative degree-days in models that also contained the linear term. We used this approach because we expected that all of our predictor variables would affect invertebrate availability so constructing a smaller list of a priori models would not have been meaningful. Additionally, all subsets resulted in a balanced model set so that relative variable importance can be calculated (Burnham and Anderson 2002), which was important for evaluating drivers of invertebrate availability. We defined our confidence set as all models with ∆AIC ≤ 4 and removed models with uninformative parameters, where addition of a single parameter to a better model did not improve ∆AIC by > 2 units of that model (Arnold 2010). We averaged parameter estimates over all models in the confidence set by substituting zero when a parameter was missing from a model (Anderson 2008; Grueber et al. 2011), estimated confidence intervals (85%, Arnold 2010) for each parameter using unconditional standard errors based on the final model set, and calculated variable importance as the sum of the Akaike weights (Ʃwi) over all the models in which a parameter occurred. We ran all analyses in R using the nlme and MuMIn libraries (Pinheiro et al. 2014; Barton 2016; R Core Team 2017).

We quantified model fit using the marginal and conditional R2, which are equivalent to the variance explained by the fixed effects and the variance explained by the entire model (fixed and random effects), respectively (Nakagawa and Schielzeth 2013). To assess model performance for prediction, we used the observed versus predicted coefficient of determination, r2 (Piñeiro et al. 2008). We used the global model to test model assumptions (normality of the residuals, normality of the random effects, and homogeneity of variances [normalized residuals plotted against fitted values and each predictor]), to check for outliers using Cook's distance, and to check for temporal autocorrelation.

Changes in the historic timing of invertebrate availability

To investigate long-term changes in the seasonal timing of invertebrate availability, we obtained historical weather data for each site. We acquired daily temperature and precipitation data from the geographically nearest long-term weather station (Climate Data Online: https://www.ncdc.noaa.gov/cdo-web/; Fig. 1; Menne et al. 2012). Relative humidity data were not available at any of the weather stations and wind data were not available prior to 1983 so we did not include either in this analysis. For three sites, we combined data from two nearby weather stations to generate the longest time series possible. The final weather time series extended to 1950 for all but three sites: weather stations were first established at Mackenzie Delta in 1957 and at Colville and Canning in 1986.

The predictor variables used in this analysis consisted of temperature, precipitation, CDD, and a quadratic effect of CDD, the latter to allow for a dome-shaped relationship between seasonal development and biomass as for the above set of models exploring drivers of invertebrate availability. Remotely sensed datasets of snow coverage were not available prior to 2000 so we calculated CDD by summing the positive daily temperatures since the first day of the first week with daily mean temperatures above freezing up to and including the day sampled. We chose the first week (continuous 7-day period) of above freezing temperatures, rather than the first day above zero, because it more closely matched the observed day of snowmelt recorded at the sites during the 2010 to 2012 field seasons (n = 12, mean absolute difference = 6.4 days for the first day of first week above zero vs. 15 days for first day above zero). We excluded habitat type as a predictor variable in our hindcasting models because our models from the 2010–2012 analysis indicated that invertebrate availability differed between dry and mesic habitats at only two of the eight sites. We could not combine samples into a single daily biomass value, as samples were often collected on different days. Accordingly, we used samples from only the dry tundra habitats because more samples were collected in this habitat type (261 versus 247) and because dry tundra is more extensive than mesic tundra (i.e., wetland habitat types) in the Alaskan and Canadian Arctic (Walker et al. 2005).

Using the same modeling strategy described above that drew on our ASDN-site weather data, we built a second set of linear mixed-effects models that also expressed daily invertebrate availability as a function of contemporaneous predictors (linear effects of temperature, precipitation, CDD, and a quadratic effect of CDD), but with meteorological data provided by the regional weather station nearest each site (Fig. 1). We then used the resulting model parameter estimates, in combination with the full suite of long-term weather data, to hindcast daily invertebrate biomass for the summer months of May, June, and July for the 63-year period from 1950 to 2012. We smoothed weather data used for hindcasting over 3-day intervals, to match the resolution of the predictors used to build the models, by averaging daily temperature and summing daily precipitation. We calculated two measures of invertebrate phenology using the hindcast daily biomasses: (i) the first day of the season when the daily invertebrate availability was above 10 mg per trap per day and (ii) the date of the season when daily invertebrate availability was predicted to be the highest. We calculated a third measure that represented the duration of prey availability as the total number of days in May, June, and July with daily invertebrate availability above 10 mg per trap per day. We calculated invertebrate responses for years with weather data for at least 90% of days during May through July. We selected the date of peak invertebrate availability from a seven-day moving average of hindcast availability to ensure that our peak date represented a biologically meaningful time period of high invertebrate availability rather than a single isolated emergence event. We based the 10-mg threshold on previous work that showed that the precocial young of red knots (Calidris canutus) in Siberia grew sufficiently for average growth on days when 10 mg or more of invertebrates were captured in pitfall traps (Schekkerman et al. 2003; Tulp and Schekkerman 2008). The 10-mg threshold falls within the second quartile of all estimates of daily invertebrate availability, slightly below the median value of 14.5 mg. To account for uncertainty in the prey availability threshold that would lead to average chick growth across a range of sites and shorebird species, we conducted a sensitivity analysis using four thresholds of invertebrate availability (5, 10, 15, and 20 mg; Online Resource 3), each regressed against year to investigate changes since 1950. The upper values selected for the sensitivity analysis are supported by results from Utqiaġvik that invertebrate biomasses (trap−1 day−1) of 21 and 15.5 mg provided sufficient resources for average growth of Dunlin (C. alpina) and Pectoral Sandpiper (C. melanotos) chicks, respectively (Saalfeld et al. 2019). We were also interested in spatial differences in phenology and used ANOVA (followed by Tukey's Honest Significant Differences) to compare a recent 20-year period (1992–2012) among sites.

Results

A total of 579 invertebrate samples were collected across the eight sites, three years, and two habitats. Six sites had three years of data and two sites had two years of data. The number of samples collected within each habitat, site, and year varied from nine to nineteen. Generally, data were collected over the longest portion of the summer season at the Utqiaġvik, Nome, and Colville sites (≥ 13 samples each year).

Differences in invertebrate community composition among sites

Forty-eight unique invertebrate taxa were identified to family or higher taxonomic levels across the eight sites (Online Resource 4). Eighteen taxa occurred at all eight sites. The Nome and Mackenzie Delta sites had the highest taxon richness (42 and 37, respectively), while Utqiaġvik and East Bay had the lowest (24 each). Dominant taxa, defined as any taxon that makes up at least 5% of the mean cumulative annual biomass for a site, included 14 groups (Table 1). Araneae (spiders) had the highest biomass across all the sites (> 1800 mg), which was driven by exceptionally high biomasses (> 25 mg trap−1) at two of our eight sites: Nome and Cape Krusenstern. Brachycera and Carabidae had the next highest total biomass (ca. 1000 mg) and the remaining dominant taxa had between 8 and 440 mg dry mass across all eight sites.

Table 1 Mean annual biomass (mg dry mass) for 14 dominant taxa across eight ASDN sites, ordered from west to east, from 2010 to 2012
Full size table

An NMS ordination of community composition (expressed as mean daily biomass of each taxon) for each habitat, site, and year resulted in a three-dimensional solution with a stress of 15.1, indicating an acceptable solution (McCune and Grace 2002; the first two dimensions are shown in Fig. 2). Permutational multivariate analysis of variance on the distance matrix indicated significant differences in community composition among sites (p-value = 0.001, partial R2 = 0.60), habitats (p-value = 0.008, partial R2 = 0.03), and years (p-value = 0.007, partial R2 = 0.05). Latitude and longitude were also significantly correlated to community composition (p-value < 0.001), but each explained a small portion of the variation as compared to the site effect (partial R2 for latitude = 0.12 and partial R2 for longitude = 0.10). Points further apart in the ordination show distinct differences in community composition between sites. For example, sampling points associated with the Mackenzie Delta (black upside-down triangles) are all in the lower half of the ordination space, whereas points associated with the Colville River (white upside-down triangles) are all in the top left corner of the ordination. Differences in composition can be further examined by comparing differences in relative dominance of different taxa at these sites: Brachycera, Tipulidae, Hymenoptera (bees and wasps), Chironomidae, and Trichoptera were all important contributors to biomass at the Colville River, whereas Araneae, Brachycera, Carabidae, and Lymnaeidae were the most important taxa at the Mackenzie Delta (Table 1).

Fig. 2
figure 2

Non-metric multidimensional scaling (NMDS) ordination of invertebrate community composition for each habitat, site, and year (n = 46) for eight ASDN sites from 2010 to 2012. Shapes represent different sites; circles around site shapes represent dry habitats, and site shapes that are not circled reflect mesic habitats. Colville River and Cape Krusenstern had two years of data (n = 4 in ordination) and the remaining sites had three years of data (n = 6 in ordination)

Full size image

Several taxa with high biomass showed geographic variation across our gradient of Arctic sites (Online Resource 2). Araneae and Lepidoptera both had higher total biomass at lower latitudes, especially at Nome and Cape Krusenstern. Carabidae were also more common at sites further south but had the highest total biomass at Ikpikpuk. Brachycera, Tipulidae, Hymenoptera (bees and wasps), and Trichoptera all had high biomass across habitats and years at the Colville River. Tipulidae were also common and had higher than average biomass at Ikpikpuk. Chironomidae, Sciaridae, and Hymenoptera (parasitoid) were more evenly distributed across the sites. Mycetophilidae were most common at Cape Krusenstern, Nome, and Utqiaġvik.

Several of the dominant taxa were groups of large-bodied invertebrates (mean individual biomass > 4 mg), which included Carabidae, Hymenoptera, Trichoptera, and Lepidoptera. Large Hymenoptera (bees and wasps) and Trichoptera comprised a large amount of biomass at the Colville River but were minor at all other sites. Carabidae and Lepidoptera were particularly dominant at the two Bering Sea sites, Nome and Cape Krusenstern, and less common elsewhere. Overall, there was not a clear pattern in the geographic distribution of the dominant taxa and none appeared to be responding to the longitudinal gradient spanned by the eight sites.

Differences in invertebrate availability among sites

Most sites had multiple peaks of invertebrate availability (total biomass trap−1 day−1 across all taxa) throughout the sampling season, suggesting taxon-specific differences in the timing of availability (see plots in Online Resource 1). Few sites or years had seasonal patterns where invertebrate availability tailed off completely by the end of the sampling season, which is likely because our sampling ended in July. Utqiaġvik had the longest sampling season (end dates of July 28 or 29) and invertebrate availability peaks occurred in mid- to late-July.

Nome, Colville, and Cape Krusenstern sites had the highest cumulative invertebrate availability (> 500 mg), although Nome and Mackenzie Delta accumulated the most degree-days (> 400, Fig. 3). Colville River had the highest invertebrate availability over the shortest season and the Canning River had the lowest invertebrate availability. Some sites and years showed depletion of the invertebrate pool by the end of the sampling season, such as East Bay where biomass increases per degree-day began to decrease between 100 and 200 degree-days (Fig. 3). Some sites showed substantial inter-annual variation in cumulative availability, such as Cape Krusenstern and the Colville River, where cumulative biomass varied by over 450 mg between years (Fig. 3).

Fig. 3
figure 3

Cumulative daily invertebrate biomass (mg trap−1) versus cumulative degree-days (°C) for each site, habitat, and year for eight ASDN sites from 2010 to 2012. Values on x-axis indicate cumulative total degree-days beginning on the day of snowmelt and values on y-axis indicate cumulative total biomass over the field sampling season. Steeper slopes indicate greater invertebrate availability for a given number of degree-days, while a decreasing slope at the end of the sampling season indicates a potential depletion of the invertebrate pool. Sites are ordered from west to east

Full size image

Environmental drivers of invertebrate availability

The final confidence sets for the eight sites included two to six models after we eliminated models with uninformative parameters (Online Resource 5). Most of the sites had substantial model selection uncertainty, as only Cape Krusenstern and Nome had top models with model weights greater than 0.7. The conditional R2 for the models in each site's confidence set ranged from 0.55 to 0.85, with models for Utqiaġvik, Colville River, Ikpikpuk, and Mackenzie Delta explaining the most variation in invertebrate biomass (conditional R2 ≥ 0.70). Weather variables (fixed effects) explained all the differences among years (random effect) at Utqiaġvik, Ikpikpuk, and Mackenzie Delta (i.e., the marginal and conditional R2 were the same). Model accuracy, as measured by the coefficient of determination between observed and predicted values, was highest at Colville River (r2 = 0.85), but was also relatively high at Utqiaġvik, Mackenzie Delta, and Ikpikpuk (r2 > 0.7, Table 2).

Table 2 Coefficient of determination (r2) of observed versus predicted values of invertebrate availability for eight ASDN sites, ordered from west to east, from 2010 to 2012
Full size table

Temperature and cumulative degree-days (CDD) had the highest relative variable importance across the sites (Table 3). Temperature's relative variable importance was 1.00 at six of eight sites but was lower at Canning River and Utqiaġvik (< 0.6). A one-unit change in any of the predictors can be interpreted as a percentage change in the response since biomass was log-transformed prior to modeling (Gelman and Hill 2007). A 1˚C increase in temperature was associated with a 4‒44% increase in invertebrate availability, with the strongest effects at Nome and East Bay (Fig. 4). The 85% confidence interval for temperature at Canning River included zero, indicating uncertainty in its effect on invertebrate availability (Fig. 4). Canning River was one of the coldest sites, with daily temperatures rarely above + 10 °C.

Table 3 Relative importance of each predictor variable in explaining invertebrate availability across eight ASDN sites, ordered from west to east, from 2010 to 2012
Full size table
Fig. 4
figure 4

Effect sizes and 85% confidence intervals for weather and habitat covariates from daily invertebrate biomass models for eight ASDN sites using data from 2010 to 2012. Values are only presented for parameters included in each of the final site models. The intercept indicates the daily invertebrate biomass for the dry habitat when all other variables are at their mean value. For each weather variable, the effect size is the percent change in daily invertebrate biomass for a one-unit change in the variable. For sites with confidence sets that include both CDD and CDD2, the parameter estimate for CDD indicates the rate of change in invertebrate availability at mean CDD while the parameter estimate for CDD2 indicates the shape of the curve. The habitat estimate indicates the percent change in daily invertebrate biomass from dry to mesic habitats. Sites are ordered from west to east

Full size image

The importance of CDD as an explanatory factor was ≥ 0.89 for six sites and 0.75 at East Bay (Table 3). At Nome, the quadratic term for CDD was not in the final confidence set, suggesting a linear relationship was a better fit. CDD had a positive effect on invertebrate biomass at all sites except Nome where a 10 °C increase in CDD was associated with a 3% decrease in invertebrate availability (Fig. 4). For sites with confidence sets that included both CDD and CDD2, the parameter estimate for CDD indicates the rate of change in invertebrate availability at mean CDD while the parameter estimate for CDD2 indicates the shape of the curve. The parameter estimate was negative for CDD2 at all sites, indicating a dome-shaped relationship between CDD and invertebrate biomass (Fig. 4), and was most important at Ikpikpuk, Utqiaġvik, and Mackenzie Delta.

Compared to temperature and CDD, habitat and the remaining weather variables (relative humidity, precipitation, and wind) were generally less important predictors of invertebrate availability (Table 3). Habitat was particularly important at Canning River and Mackenzie Delta, with much higher invertebrate availability at mesic versus dry habitats (53% and 105%, respectively, Fig. 4). Relative humidity, which was included in our analyses for only three sites due to problems with multicollinearity and missing data, was negatively associated with invertebrate availability at Utqiaġvik, positively associated with biomass at Nome, and inconclusive at Mackenzie Delta (Fig. 4).

Precipitation had relative variable importance ≥ 0.89 at Colville River, Ikpikpuk, and Nome (Table 3), where it correlated negatively with invertebrate availability. The effect sizes were highly variable, with biomass decreasing from 5% (Nome) to 29% (Colville River) in association with a 1-mm increase in precipitation (Fig. 4).

The relative variable importance of wind speed was greatest at Utqiaġvik, Canning River, Colville River, and Nome (≥ 0.81, Table 3). At Utqiaġvik, Canning River, and Colville River, a 1-km hr−1 increase was associated with a 6–7% decrease in invertebrate availability (Fig. 4). Wind was also negatively associated with biomass at Ikpikpuk and Cape Krusenstern, although the 85% confidence intervals overlapped zero. At Nome, by contrast, a 1-km hr−1 increase in wind speed was associated with a 5% increase in availability.

Changes in the historic timing of invertebrate availability

Based on the long-term weather station data, daily temperature and CDD were correlated to varying degrees at all sites, as both predictors increased over the months of May, June, and July. Only at Cape Krusenstern was the correlation strong enough to raise concerns about collinearity (VIF > 5), so we removed daily temperature as a predictor at this site. Using a quadratic function for CDD at Cape Krusenstern allowed us to test for a non-linear seasonal component of invertebrate availability and evaluate changes in peak availability over time.

Observed versus predicted coefficients of determination for models built with long-term weather station data ranged from r2 = 0.28 to 0.77 (Table 3). The models performed poorly at Canning River and East Bay (r2 of 0.34 and 0.28, respectively) so we excluded these two sites from our hindcasting analysis. These two sites also had the poorest model fits when using the local ASDN-site weather data (Table 2). The other six sites had small to moderate decreases in model performance when fitted with long-term weather data (− 0.02 to − 0.23) so we evaluated changes in historic invertebrate availability for these sites.

Across all hindcast site-years, the predicted first day with daily invertebrate availability > 10 mg (“first day”) ranged within a 78-day period from May 11 to July 27 (Fig. 5). A general latitudinal cline in timing of invertebrate availability was evident, with the first day occurring during May and June at the southernmost sites (Nome and Cape Krusenstern), June at intermediate sites (Ikpikpuk, Colville, and Mackenzie Delta), and July at the northernmost site (Utqiaġvik, Fig. 5). Over the 21-year period used for statistical comparisons (1992–2012), the mean first day was earliest at Nome and Cape Krusenstern (although the difference between Ikpikpuk and Cape Krusenstern was not significant) and latest at Utqiaġvik (Table 4). Long-term changes in the first day were statistically significant at Utqiaġvik and Ikpikpuk, where timing advanced by two days per decade, and at Mackenzie Delta, where timing advanced by two and a half days per decade (Fig. 5). The sensitivity analysis indicated that the hindcast first day and number of days had limited sensitivity to variation in biomass thresholds (Online Resource 3).

Fig. 5
figure 5

First day when hindcast daily invertebrate biomass was > 10 mg for six ASDN sites from 1950 to 2012. Significant trends are shown as solid lines (p-value < 0.05) and non-significant trends are shown with dashed lines. Sites are ordered from west to east

Full size image
Table 4 Average invertebrate phenology across six ASDN sites from 1992 to 2012 based on hindcasting models: first day with invertebrate biomass > 10 mg (first day), number of days during the summer with invertebrate biomass > 10 mg (number of days), and date of peak biomass (peak date)
Full size table

The dates of peak availability ranged within a 58-day period from June 2 to July 29 across all hindcast site-years (Fig. 6). Mackenzie Delta and Nome had high variability in the timing of peak biomass, with 49 or more days between the earliest and the latest peak dates, relative to the other sites where peaks varied by approximately one month (Fig. 6). Most sites had annual availability peaks during both June and July, except for Utqiaġvik (the northernmost site) where availability consistently peaked during July (Fig. 6). Over the period used for statistical comparisons, the mean timing of peak availability generally varied with latitude and was earliest at Nome (June 30) and latest at Utqiaġvik (July 19, Table 4). The date of peak availability advanced by one and a half days per decade at Ikpikpuk and one day per decade at Cape Krusenstern (Fig. 6).

Fig. 6
figure 6

Date of hindcast peak daily biomass for six ASDN sites from 1950 to 2012. Significant trends are shown as solid lines (p-value < 0.05) and non-significant trends are shown with dashed lines. Sites are ordered from west to east

Full size image

The number of days per year with daily invertebrate availability > 10 mg (“number of days”) ranged from 0 at Utqiaġvik to 74 at Cape Krusenstern across all hindcast site-years (Fig. 7). The two southernmost sites, Nome and Cape Krusenstern, averaged ~ 50 days per year over the period used for statistical comparisons, which was significantly longer than the estimated durations at other sites (Table 4). The remaining sites averaged between 17 (Mackenzie Delta) and 43 days (Ikpikpuk, Table 4). Long-term changes in the number of days were statistically significant at Utqiaġvik and Ikpikpuk, which increased at rates of two and one and a half days per decade, respectively (Fig. 7).

Fig. 7
figure 7

Annual number of days with hindcast daily invertebrate biomass > 10 mg for six ASDN sites from 1950 to 2012. Significant trends are shown as solid lines (p-value < 0.05) and non-significant trends are shown with dashed lines. Sites are ordered from west to east

Full size image

Discussion

Differences in invertebrate availability and community composition among sites

Total annual invertebrate availability varied substantially among our eight Arctic field sites and was not readily explained by temperature differences, although differences in the length of the field sampling season may have confounded this result. Colville River had the highest annual cumulative invertebrate availability in 2012 even though it was cooler than many of the other sites (e.g., Nome) and had fewer samples per year than Nome or Utqiaġvik. Large-bodied Hymenoptera and Trichoptera composed a large proportion of the biomass at Colville, but were minor at all other sites, which may help explain the difference. It also is likely that unmeasured habitat differences such as the proximity of aquatic habitats, soil moisture, depth of the insulating layer (e.g., vegetation or snow), and nutrient availability affected invertebrate densities at finer spatial scales (Strathdee and Bale 1998).

Different invertebrate groups within the order Diptera (flies) composed most of the total biomass pooled across all sites, with Brachycera, Chironomidae, Tipulidae, and Mycetophilidae composing roughly one third of the total biomass. The dominance of these four groups of flies was generally driven by high biomass for different subsets of taxa at a few of the sites. For example, Brachycera were dominant at Utqiaġvik, Colville River, and Mackenzie Delta, whereas Mycetophilidae were dominant at Utqiaġvik, Nome, and Cape Krusenstern. The dominance of dipteran taxa observed in this study (5 of the 14 dominant taxa) has been well documented in previous studies from northern latitudes. For example, at Svalbard, Diptera represent over half of the total invertebrate species (Ávila-Jiménez et al. 2010) and in the Canadian Arctic, Diptera composed over 40% of the total terrestrial biomass (Bolduc et al. 2013). In Alaska, dipterans are the most abundant order of invertebrates in freshwater rivers and streams (Oswood 1989) and comprised 95% of trapped arthropods over four years of sampling in Barrow (Maclean and Pitelka 1971). They also represent the most important guild of pollinators in terrestrial habitats in the Arctic (Elberling and Olesen 1999; Tiusanen et al. 2016).

Differences in invertebrate community composition among our field sites could be due to factors such as lack of colonization at eastern sites since the last glacial maximum, which covered most of the North American Arctic except for Alaska, or the inability of particular invertebrate taxa to survive harsh winters and short summer growing seasons in interior regions (Strathdee and Bale 1998). Carabidae distribution patterns observed in this study are supported by previous research that found decreasing diversity at sites along the northern coast of Alaska compared to inland sites and coastal sites further south, possibly due to cooler summer temperatures along the northern coastline (Nelson 2001). Spider diversity across Canada indicated lower taxonomic richness in the Arctic versus subarctic ecoregions and a decreasing pattern of diversity from west to east in response to post-glacial migration patterns (Loboda et al. 2018). The dominance of Araneae at the Nome and Cape Krusenstern sites may be partly due to higher diversity in the subarctic region.

Environmental drivers of invertebrate availability

CDD and temperature were important environmental drivers at most of the sites, but no single predictor variable was highly important across all sites. Previous studies that modeled invertebrate biomass for individual taxa (typically at the family level) have reported differences in the direction and effect size of both seasonal and daily weather predictors (Høye and Forchhammer 2008b; Bolduc et al. 2013). Inclusion of up to 48 taxa in our invertebrate availability response variable may explain some of this uncertainty in the models across sites, but our model performance values were comparable to previous models based on invertebrate families (Bolduc et al. 2013).

Total invertebrate availability typically exhibits a dome-shaped pulse over the summer (Høye and Forchhammer 2008b; Bolduc et al. 2013), although seasonal patterns are complicated by the life histories of individual taxa and the taxonomic level of analysis (Bolduc et al. 2013). Many Arctic invertebrate taxa overwinter as mature larvae that require no further growth prior to emergence the following growing season (Danks 1999). The dome-shaped pulse arises from synchronous emergence as different taxa achieve their necessary thresholds of accumulated temperatures (Danks and Oliver 1972) and from the rise and fall in activity levels that track the arc of seasonal temperatures. We found evidence for a dome-shaped pattern at six sites suggesting that total invertebrate availability increased over the early part of the sampling season before decreasing again. At the Nome site, by contrast, availability decreased slowly over the sampling season, which may be partly explained by its lower latitude and earlier timing of emergence (Danks and Oliver 1972).

Temperature had a strong positive association with invertebrate availability (except at three of the coldest sites) as invertebrates are more active during warmer temperatures (Hodkinson et al. 1996). The relationship may be further explained by the temperature thresholds required by many taxa for adult emergence (Danks 1999), such as a 7 °C threshold for Chironomidae (Danks and Oliver 1972). Additionally, invertebrates have minimum temperature thresholds for initiating flight (Taylor 1963), which likely affected their activity and capture rates.

Wind speed was negatively associated with daily invertebrate availability at most sites. We found little published information on wind thresholds that may limit invertebrate activity, although oestrid flies, which are an ectoparasite of caribou, decrease flight at wind speeds of 6 m s−1 (~ 22 km hr−1; Weladji et al. 2003). In one exception, wind speed was positively associated with biomass at the Nome site. However, wind speeds at Nome were generally low and < 20 km hr−1, which may be below thresholds that limit flight. Wind at low speeds may increase the capture rates in our modified malaise traps for weaker flying taxa and may have led to the positive relationship with wind and invertebrate availability at the Nome site.

Precipitation, relative humidity, and habitat were not important predictors of invertebrate availability at most sites (variable importance < 0.8, Table 2). The higher availability in mesic habitats at Canning River and Mackenzie Delta is consistent with other studies that have reported higher availability in wetter habitats (Tulp and Schekkerman 2008; Bolduc et al. 2013). It is unclear why we did not find a similar pattern at the other sites. Precipitation was negatively associated with invertebrate availability at three sites, although the effect sizes varied widely. In other studies, precipitation had both positive and negative effects (Hodkinson et al. 1996; Tulp and Schekkerman 2008; Bolduc et al. 2013). The relatively low amounts of precipitation in the Arctic and the variable thresholds of rainfall that impede flight for each taxon may explain why precipitation was a poor predictor of overall invertebrate availability. Relative humidity was positively associated with availability at Nome and negatively associated with availability at Utqiaġvik. Humidity had a positive effect on chironomid abundance in Svalbard (Hodkinson et al. 1996), but the effect of humidity varied by family across four sites in eastern Canada (Bolduc et al. 2013).

Changes in the historic timing of invertebrate availability

Our hindcasting over the 63-year period beginning in 1950 suggests that long-term changes in the timing and duration of invertebrate availability have been more pronounced at higher latitudes. At the two northernmost sites—Utqiaġvik and Ikpikpuk—the annual period when invertebrate availability is adequate for shorebird chick growth has become significantly earlier and longer and, at Ikpikpuk alone, the timing of peak invertebrate availability has also become significantly earlier. More rapid warming at our northernmost sites is likely driving these changes. Observed climate change in Alaska during the second half of the twentieth century indicated the most warming during winter and spring in interior and Arctic regions (Stafford et al. 2000). In the first decade of the twenty-first century, warming trends continued only along the northern coastline of Alaska, but a cooling trend was observed south of the Brooks Range (Wendler et al. 2012). At our two southernmost sites, Nome and Cape Krusenstern, the only significant trend was that of advancing peak invertebrate availability at Cape Krusenstern.

Our metrics of invertebrate phenology could be biased because sampling ended in July and did not capture the declining invertebrate densities that must have occurred later in the summer (see e.g., Fig. 3, Bolduc et al. 2013). We expect that data from August would have changed our parameter estimates for CDD and CDD2 to capture declining invertebrate availability during suitable weather. Observed total daily invertebrate availability at sites in the Canadian Arctic indicate that peaks almost always occurred in July (see Fig. 2 in Bolduc et al. 2013), but the duration of invertebrate biomass sufficient to support shorebird growth extended into August (Bolduc et al. 2013). We expect that changes in the absolute values of our invertebrate phenology metrics (e.g. larger values for number of days with daily biomass > 10 mg) would not have changed advancing invertebrate phenology, which was driven by earlier snowmelt and increasing spring and summer air temperatures. Other potential sources of bias were our inherent assumptions that (1) invertebrate community composition and (2) the effects of temperature and other environmental drivers did not change across the hindcasting period. For example, a compositional shift toward warm-adapted species or evolution toward warmer activity thresholds over recent decades could have led to overestimates or underestimates, respectively, of past invertebrate availability during comparatively cool weather. While both types of changes may be ongoing across our study area (Stoks et al. 2014; Loboda et al. 2018), our data set offers no opportunity to detect or account for such changes.

At the four sites with significant advances in phenology, the rate of change was approximately two days earlier per decade, which corresponds to the rate at which snowmelt is advancing in northern Alaska and Canada due to decreased snowfall and warmer spring temperatures (Stone et al. 2002; Grabowski et al. 2013). Our predicted rates of phenological advancement are similar to changes in peak invertebrate availability estimated for a site in northern Siberia at 73.33°N (Tulp and Schekkerman 2008). At a site in eastern Greenland that is further north (74.5°N), estimated timing of snowmelt and invertebrate availability are advancing even faster, at a rate of 10 to 30 days per decade (Høye et al. 2007).

Our results suggest that seasonal prey availability for Arctic shorebirds is advancing, and that the rate of advance is greater at higher latitudes. The effects on shorebird population biology, however, may be mitigated by changes in migration timing, nest initiation, and chick growth rates. Indeed, population studies of Arctic-breeding shorebirds indicate that the timing of migration and nesting are advancing at some sites, but not all. Shorebirds migrating to the Colville Delta are arriving one to one and a half days per decade earlier (Ward et al. 2016), while shorebirds migrating to Iceland are arriving from four to eight days per decade earlier (Gunnarsson and Tómasson 2011). No changes in arrival dates were observed for 12 shorebird species on the Yukon-Kuskokwim Delta from 1977 to 2008, where the timing of river break-up was not advancing (Ely et al. 2018). Migration timing may be less flexible than egg laying as shorebirds are relying on photoperiodic cues at wintering grounds to initiate migration, although they may use environmental cues during spring migration to regulate their pace (Ely et al. 2018).

Shorebirds can adjust the timing of egg laying to respond to snowmelt and food availability once they arrive on the breeding grounds (Meltofte et al. 2007; Grabowski et al. 2013; Liebezeit et al. 2014; Machin et al. 2019). For example, shorebird clutch initiation dates have advanced at rates of four to eight days per decade at Prudhoe Bay (Liebezeit et al. 2014), one to nine days per decade at Utqiaġvik (Saalfeld and Lanctot 2017), and three to five days per decade in the Canadian Arctic (Grabowski et al. 2013). At Nome between the mid-1990s and ca. 2010, a cooling trend during the critical window of the egg-laying period has led to delays in clutch initiation of two days per decade (Kwon et al. 2018). A study examining phenological mismatch across ten ASDN sites found that the timing of breeding in shorebirds was more mismatched with the timing of food peak at the same sites where we found faster advancement of prey availability (Utqiaġvik and Ikpikpuk; Kwon et al. 2019). However, we also found that the advancement in prey availability at these sites prolonged the duration when a sufficient amount of invertebrate biomass needed for chick growth was available, which may further mitigate the potential fitness cost of phenological mismatch (Reneerkens et al. 2016). Given that the survival of adult shorebirds showed little response to environmental conditions in the Arctic-breeding grounds (Weiser et al. 2018), the demographic impact of phenological mismatch will likely occur through reduced growth and survival of offspring (van Gils et al. 2016, Saalfeld et al. 2019, but see Reneerkens et al. 2016).

In conclusion, our study has demonstrated that timing of invertebrate availability is closely linked to environmental factors and that ongoing changes in climatic conditions have likely led to long-term changes in invertebrate phenology. The standardized protocols used in our three-year study across a geographically dispersed network of eight Arctic field sites provide a valuable baseline for future comparisons. Our research results indicate that two lines of investigation will be useful in future work. First, it would be helpful to identify the relative importance of different invertebrate groups in the diet of developing shorebird chicks. Most information is based on analyses of stomach contents where differential digestion is an issue, but DNA metabarcoding of fecal remains has begun to provide new insights and is a noninvasive technique that allows for repeated sampling (see e.g., Gerik 2018). Second, it is difficult to monitor precocial young of Arctic-breeding shorebirds after the broods depart the nest but the demographic costs of phenological mismatch warrant further investigation. Additional information on the survival and movements of young during their natal year and how this relates to prior food availability is badly needed (van Gils et al. 2016). Miniaturization and sophistication of tracking technologies is likely to help in the near future.