Willow drives changes in arthropod communities of northwestern Alaska: ecological implications of shrub expansion

. Arthropods serve as complex linkages between plants and higher-level predators in Arctic ecosystems and provide key ecosystem services such as pollination and nutrient cycling. Arctic plant communities are changing as tall woody shrubs expand onto tundra, but potential effects on arthropod abundance and food web structure remain unclear. Changes in vegetation structure can alter the physical habitat, thermal environment, and food available to arthropods, thereby having the potential to induce cascading effects throughout the ecosystem. We evaluated relationships between the abundance, biomass, and community composition of arthropods and the cover of several shrub taxa across tundra – shrub gradi-ents in northwestern Alaska. While previous research had found a general positive association between arthropod biomass and shrub cover, we found heterogeneity in this relationship with ﬁ ner-scale examina-tion of (1) shrub taxa, (2) arthropod taxa, and (3) arthropod guilds. Abundance and biomass of arthropods showed strong, positive associations with the amount of cover of willow ( Salix spp.) but were not signi ﬁ cantly in ﬂ uenced by shrub birch ( Betula spp.) or ericaceous (Ericaceae) vegetation. Signi ﬁ cant shifts in arthropod community composition were also associated with willows. Among trophic groups of arthropods, herbivores and pollinators were most positively associated with willow cover. Due to geographical variation in both dominant shrub taxa and their rates of expansion, effects on arthropod communities are likely to be heterogeneous across the Arctic. Taken together, our results suggest that shrub expansion could increase food availability for higher-level insectivores and shift Arctic food web structure.


INTRODUCTION
Shrub cover is increasing in the Arctic (Tape et al. 2006), potentially modifying animal community structure and therefore ecosystem function (Christie et al. 2015, Tape et al. 2016. Warmer, longer summers create conditions favorable for the colonization of tundra habitats by woody shrubs, particularly willow (Salix spp.), birch (Betula spp.), and alder (Alnus spp.) (Racine et al. 2004, Tape et al. 2006. Shrub thickets have higher biomass but lower plant diversity than graminoid tundra and differ in abiotic factors such as solar radiation, soil moisture, soil temperature, and snow cover (Myers-Smith et al. 2011, Lantz et al. 2012, Sweet et al. 2014). These biotic and abiotic factors may be contributing to recent changes in the distribution and abundance of Arctic and subarctic animals (Rich et al. 2013, Christie et al. 2015, Tape et al. 2016. Understanding the impact of shrub expansion on Arctic arthropods is particularly important because arthropods comprise much of Arctic biodiversity and provide many ecosystem services, including pollination, nutrient cycling, decomposition, predation, and food for insectivores (Hodkinson et al. 2013). Arthropods are of particular importance to many migratory birds that depend on them as a protein source for nestlings (Schekkerman et al. 2003. Predicting the impacts of shrub expansion on Arctic arthropods remains difficult, however, because increased shrub dominance may exert counteracting forces on arthropods in cold, northern climates. Shrub thickets trap snow, which provides insulation and may increase winter survival for many arthropod taxa, but shrubs also reduce solar radiation, which delays snowmelt and lowers spring soil temperatures (Lantz et al. 2012, Sweet et al. 2014, potentially delaying development (Gilbert and Raworth 2000). Shrub thickets provide more vertical habitat, which may promote arthropod abundance and taxonomic richness (Koricheva et al. 2000, Rich et al. 2013; however, shrub expansion is associated with decreased plant diversity (Walker et al. 2006), which in turn could decrease arthropod diversity (Siemann et al. 1998). To better predict the ecological effects of shrub expansion on Arctic ecosystems, we need more empirical data on how arthropod communities differ between tundra and expanding shrub habitats.
While prior research in Arctic Alaska has found higher arthropod abundance and biomass in shrub thickets compared with open tundra (Rich et al. 2013 and uncovered some differences in arthropod composition at the order level (Rich et al. 2013), these relationships are likely to vary among arthropod guilds and lower-level taxa. For example, the greater plant biomass available in shrub thickets may result in higher abundance of generalist herbivores but not specialists whose hosts are excluded from shrubs (Strong et al. 1984). Some pollinators may be negatively associated with shrubs because many species of flowering plants are found only in open tundra (Swanson et al. 1985). In addition, shrub taxa vary in the resources they provide to arthropod herbivores, pollinators, and associated natural enemies (Strong et al. 1984, MacLean and Jensen 1985, Mulder 1999. While grouping plants by growth form (e.g., deciduous shrubs vs. graminoids) is a common and useful practice in studies of the ecological impacts of shrub expansion (Elmendorf et al. 2012, Rich et al. 2013, it leaves open questions about important biological differences among shrubs, such as height and palatability (Christie et al. 2015, Saccone et al. 2017. Given differences among shrub taxa in the rate and spatial extent of expansion (Myers-Smith et al. 2011, Elmendorf et al. 2012, it is important to understand how arthropod community composition varies with shrub taxa. We examined patterns of arthropod abundance, biomass, and community composition relative to the structure and taxonomic composition of vegetation as part of a study of insectivorous, migratory birds in the tundra−shrub transition zone of northwestern Alaska (Thompson et al. 2016). We classified arthropods to a low enough taxonomic level to infer the trophic function of most specimens. Based on samples from a series of plots at coastal and inland sites, we constructed a set of models to determine which shrub types most strongly influenced the arthropod abundance, biomass, and composition at a given site. We also analyzed slope, representing sun and wind exposure, since there is evidence that this variable helps explain variation in arthropod abundance and diversity in the Arctic (Hansen et al. 2016). Based on previous studies in Arctic Alaska (Rich et al. 2013, we expected that numerical abundance and biomass of arthropods would generally be higher in areas with greater shrub cover and height. We were particularly interested in extending previous findings by determining whether arthropod abundance, biomass, or composition would also vary depending upon the taxonomic identity of the dominant shrubs and whether responses were specific to arthropod taxonomic groups or feeding guilds. We expected that palatable shrub species with few defense compounds (e.g., willow) would have a more positive association with arthropod abundance and biomass than less palatable shrub species (e.g., ericaceous shrubs). Finally, v www.esajournals.org we expected that different taxonomic groups and guilds of arthropods would vary in their response to overall cover and taxonomic identity of dominant shrubs, with the biomass of herbivorous arthropods likely to be positively associated with shrub cover due to increased biomass of a food source. Our study contributes valuable empirical data on arthropod communities in a rapidly changing Arctic and provides a more nuanced understanding of how shrub types vary in their effects on arthropod assemblages.

Study area
We conducted our study on the Seward Peninsula in northwestern Alaska, USA (Fig. 1). The peninsula encompasses the transition zone between boreal forest and tundra and is currently undergoing rapid climate-mediated shrub expansion into former tundra vegetation (Silapaswan et al. 2001, Racine et al. 2004. Between 1986 and 1999 in this region, shrubs advanced about 100 m to the northwest and increased in cover (Silapaswan et al. 2001).
We established a series of 600 × 600 m plots at two study sites, one interior and one coastal (Fig. 1). In 2015, we established 12 plots in pairs (<4 km apart) that contrasted in dominance of tundra vegetation (Fig. 1a) vs. taller shrubs (Fig. 1b), with three pairs of plots at each study site. In 2016, we replaced one pair of plots at the interior study site with a new single tundra-dominated plot due to logistical constraints. Most of the study plots had heterogeneous vegetation, with shrubs and tundra interspersed; therefore, all of our analyses consider vegetation characteristics at finer spatial scales than whole plots.

Arthropod sampling
In 2015 and 2016, we collected arthropods weekly during the warmest period of the year when arthropod abundance and activity are highest (8 June-28 July 2015; 3 June-27 July 2016). We used pitfall trapping (targeting ground-dwelling arthropods) and sweep-netting (targeting canopy-dwelling arthropods) to sample the above-ground arthropod community in each plot. We established three pitfall arrays 200 m apart in a line through the center of each plot, where each array was considered a sampling location. Arrays consisted of three cups (7.6 cm diameter; 16.5 cm tall) placed in a line 1 m apart, buried flush with the soil surface, filled with 2 cm of propylene glycol, and covered with a plywood square 3-4 cm above the rim to protect the cup from rainwater. Contents of the three cups (subsamples) in each array were combined into a single sample for identification and enumeration. In a few cases when pitfalls had been open for 6 or 8 d instead of 7, the counts of insects caught were normalized to 7 trap-days. There were seven weekly sampling periods each season.
Sweep-net transects consisted of 50 side-toside sweeps over 25 m in a direction randomly chosen for each transect throughout the sampling period from the same sampling locations used for the pitfall arrays (i.e., central cup of each array). Sweeps were aimed at the top 25% of vegetation within a 2-m swath along each transect. For most plots, we collected sweep-net samples at two (in 2015) or all three of the pitfall sampling locations (in 2016) per plot each week. In three plots with highly heterogeneous vegetation, we conducted sweep-net samples at five collection points in 2016 (the three pitfall sampling locations and two additional locations 200 m from the central sampling location).
Using published keys (Triplehorn andJohnson 2005, Marshall 2006), we sorted arthropod samples and identified specimens to family, although some were grouped at higher taxonomic levels due to similarity in trophic function (Lepidoptera, Araneae) or the difficulty in distinguishing families (Diptera: Muscomorpha). While lower classifications would be preferable, this approach helped to ensure accuracy in identifications while still allowing us to analyze functional differences among arthropod groups. In addition, published keys are not available for some families of Arctic arthropods, making classifications below family difficult if not impossible in some cases. We identified the trophic guild of each taxon (Triplehorn andJohnson 2005, Marshall 2006) and when species within families varied in their dietary habits we categorized the trophic guild as 'varies within family' (Appendix S1: Table S1).
We aggregated counts of arthropods across the season for each sampling location, so each data point represented the sum of 7 (pitfall) or 8 v www.esajournals.org (sweep) samples taken periodically at one sampling location. Pitfall trapping and sweep-netting are composite measures of abundance and activity, and while vegetation structure can affect arthropod activity, aggregate counts across a season can provide more reliable abundance estimates for some taxa (Baars 1979). To estimate arthropod biomass, we dried as many individuals as possible (Appendix S1: Table S2) of each taxon of arthropods for 48 h at 40°C, calculated mean individual dry mass, and multiplied this by counts of individuals to estimate biomass of each sample. When size was highly variable within a taxon, we separated individuals by length and obtained mean dry mass for each size class. Individuals within each taxon were counted and weighed separately by life stage (larva, nymph, adult), but life stages were combined for analyses of community composition.

Vegetation sampling
We sampled vegetation during a 3-week period in mid to late July 2015 (initial plots) or 2016 (added plot), after most tundra plants had reached peak greenness ; we assumed changes in vegetation cover would be minimal between consecutive summers. This study was part of a larger initiative focused on climate effects on Arctic migratory birds and vegetation, and therefore vegetation measurements represent a balance between characteristics and spatial scales relevant to birds and arthropods. We characterized vegetation within five 0.5 × 0.5 m quadrats, 5 m apart along one 20-m, randomly oriented transect that was centered on each arthropod sampling location (where pitfall arrays were located and all sweepnet transects started). Plant diversity is generally low in our study area, and this spatial scale provides a good representation of the heterogeneity of the surrounding plant community relevant to ground-dwelling and canopy-dwelling arthropods, although flying insects in particular are likely using habitats beyond the spatial extent of our vegetation sampling. We used modified Daubenmire (1959) cover class categories (0%, 1-5%, 6-25%, 26-50%, 51-75%, 76-95%, and 96-100%) to characterize vegetation cover. Within each quadrat, we estimated cover classes for the following shrub types: willow (Salix alaxensis, S. glauca, S. pulchra, and S. richardsonii), birch (Betula nana and B. glandulosa), alder (Alnus crispa), and ericaceous shrubs (primarily Vaccinium uliginosum, V. vitis-idea, Empetrum nigrum, Cassiope tetragona, Ledum decumbens, and Arctostaphylos alpina). The four willow species of interest are not the only willows present on our sites, but they are the only willow species that have the potential to grow into tall thickets (over 30 cm) and therefore are the most relevant to studying the effects of shrub expansion. We grouped ericaceous taxa because all are slowgrowing and short (usually <30 cm), are primarily associated with tundra habitats, and do not play a role in ongoing shrub expansion. We measured height of each shrub type as the mode in cm of 5-10 haphazardly dropped points within each quadrat. To determine percent vegetation cover in a quadrat, we converted cover classes to the midpoint percentage (0%, 2.5%, 15%, 37.5%, 62.5%, 85%, and 97.5%). Shrub cover and height were averaged across quadrats within a transect (5 per transect) for the overall transect value used in analyses. We excluded alder from analyses because it was observed on only one transect. Slope, a potentially important abiotic variable related to sun and wind exposure, was measured for each point with a clinometer over a horizontal distance of 10 m. Elevation was fairly consistent among plots (range: 42-287 m) and was not included in our analyses of arthropod abundance and biomass.
All data have been archived and are publicly available in McDermott et al. (2021). Arthropod voucher specimens are in the University of Alaska Museum of the North (accession number UAM-2014.24-USGS-Ento).

Statistical analysis
We built on previous research that showed a positive relationship between overall arthropod abundance and shrub dominance in the Arctic (Rich et al. 2013 by examining possible heterogeneity of this effect among shrub taxa. Toward this end, we chose to examine the impact of major shrub groups present in our study region (willow, birch, and ericaceous shrubs-again, alder was not analyzed since it was observed on only one transect) as well as slope. To facilitate comparison between models examining different taxonomic and trophic groups of arthropods, we used the same set of predictor variables for all models of arthropod abundance and biomass. We conducted all statistical analyses using R version 3.5.1 (R Development Core Team 2018). We built linear mixedeffects models and generalized linear mixed models using package lme4 (Bates et al. 2015). Sweep-net samples (n = 63 points) and pitfall samples (n = 69 points) were analyzed separately. To account for repeated measures across years and the unbalanced structure of our v www.esajournals.org experimental design, we included year and site as fixed effects (since each has only two levels; Bolker et al. 2009) and plot as a random effect. We examined community composition using non-metric multidimensional scaling (NMDS) with the R package vegan (Oksanen et al. 2018).
Data preparation.-We analyzed shrub variables and arthropod data at the level of sampling locations within a plot. To address potential multicollinearity between shrub variables, we calculated Pearson correlation coefficients for the suite of candidate predictors in a model. Height and cover for each shrub type were highly collinear. In each linear model, we decided to model cover as our predictor variable since it is the shrub attribute measured most commonly in studies of shrub expansion, many of which are done via remote sensing. In addition, unlike for height, values of zero for cover are biologically meaningful. With height excluded, none of our remaining predictors had correlation coefficients >0.4, with the exception of linear and quadratic terms for the same variable. We centered and scaled all predictors to aid in model convergence.
Overall abundance and biomass. -Using the model structure described above, we analyzed the influence of shrub characteristics on overall arthropod abundance and biomass while accounting for abiotic conditions. If the dispersion parameter (ϴ) of any negative binomial model was close to 1, we compared the fit of a Poisson model and selected the model with the lower Akaike's information criterion (AIC) score. For all arthropod taxa, models of overall abundance were best fit by a negative binomial distribution, compared with a Poisson distribution or zero-inflated models, as determined by comparing AIC scores. Biomass data were normally distributed after transformation following the Box-Cox procedure (pitfall samples were cuberoot transformed, sweep-net samples were logtransformed; Shapiro-Wilk test: pitfall, P = 0.17; sweep-net, P = 0.51).
To allow for potential non-linear relationships of arthropods with shrub taxa, we added a quadratic term to the model for each potentially nonlinear variable. Quadratic terms were kept if they improved model fit as assessed by AIC. To address the influence of outliers, we plotted the residuals against fitted values and each predictor variable. Observations with residuals of SD >3 were checked to ensure accuracy, and if they represented real data points, they were truncated to the next highest value and models re-fit. We opted to retain and truncate, rather than omit, these values (following Thompson et al. 2016) because they represented real data points of dense clusters of patchily distributed arthropods that, if unadjusted, could bias model estimates. Model assumptions were again checked. Standard goodness-of-fit measures for GLMM models are lacking, as random effects make it unclear how to define what proportion of variance is explained by the model (Bolker et al. 2009). To provide some measure of how well our models fit the data, we calculated marginal and conditional pseudo-R 2 following Nakagawa and Schielzeth (2013). In addition, we fit a null model with random effects only (including year and site, which while coded as fixed effects are conceptually analogous to random effects). We report ΔAIC from null to full model.
Taxon-specific abundance. -To understand the relationship between taxonomic abundance and shrub cover, we built models for the six most abundant arthropod orders (Araneae, Coleoptera, Diptera, Hemiptera, Hymenoptera, and Lepidoptera; together representing 97% of total abundance). We analyzed counts from only the sample type (pitfall, sweep-net) in which most individuals (over 60%) were collected. Hymenoptera were abundant (more than 40%) in both pitfall traps and sweep-nets, so we performed separate analyses for each sample type. Most models were best fit by a negative binomial distribution, with the exception of Lepidoptera, for which we used a Poisson distribution. Model-fitting and checking proceeded as outlined above.
Biomass of trophic levels.-We were interested in how shrub dominance affects trophic structure, and for this analysis we chose to analyze arthropod biomass as this measure may be more important to the flow of energy across trophic levels than numerical abundance. To understand the relationship between arthropod biomass and shrub cover, we built separate models for herbivores, predators, and pollinators (the three most common trophic levels). For some very rare taxa, we did not have biomass data and therefore these taxa were not included in this analysis. In cases where a taxon served multiple trophic functions (for example, Cantharidae are both v www.esajournals.org pollinators and predators), we included them in both trophic groups for modeling. Herbivores comprised Chrysomelidae, Curculionidae, Elateridae, Silphidae, Cecidomyiidae, Aphidoidea, Cicadellidae, Delphacidae, Lygaeoidea, Miridae, Psyllidae, Lepidoptera larvae, and Tenthredinidae larvae. Predators comprised Araneae, Cantharidae, Carabidae, Staphylinidae, Dolichopodidae, Empididae, Rhagionidae, Tenthredinidae adults, and Neuroptera. Pollinators comprised Cantharidae, Bibionidae, Ceratopogonidae, Chironomidae, Culicidae, Simuliidae, Empididae, Syrphidae, Lepidoptera adults, and Trichoptera (percent contribution of each taxonomic group to overall biomass of each guild is listed in Appendix S1: Table S3). As above, we analyzed biomass from only the sample type (pitfall or sweep-net) that represented the majority of biomass (more than 60%) collected for each trophic level. Pitfall traps caught the majority of biomass of predators, whereas sweep-net sampling caught the majority of biomass of herbivores and pollinators. Some taxa were only present in the other sampling type (e.g., 100% of Apidae, which are pollinators, were caught in pitfall traps) and are therefore not listed in the categories above nor were they analyzed, but all such groups represent a small fraction (<1%) of overall biomass. Appropriate transformations were identified using the Box-Cox procedure. Model-fitting and checking proceeded as outlined above.
Community composition.-To examine how vegetation characteristics influenced the arthropod community as a whole, we used NMDS with the Bray-Curtis index of dissimilarity to visualize community similarity. Issues with convergence on a two-dimensional solution were resolved by setting k = 3. We conducted permutation tests to determine the relative significance of shrub cover and height on community composition, with plot as a blocking factor to account for spatial similarity.
Willow species were present on most transects and were typically <50 cm tall, but a few transects included very tall thickets (>3 m). Ericaceous plants were present on almost every transect but tended to be short (<30 cm). Birch shrubs were present on just over half of transects, usually <20 cm tall but occasionally up to 1 m in height.

Overall arthropod abundance
Overall abundance of arthropods increased significantly with willow cover for both pitfall and sweep-net samples ( Table 2, Fig. 2a, b). Neither shrub birch nor ericaceous cover was significantly associated with overall abundance in either sample type (Table 2, Fig. 2c-f). Arthropod abundance in pitfall traps was negatively associated with slope, but arthropod abundance in sweep-net samples had no significant relationship with slope. Abundance of pitfall taxa was greater in 2016 and in the coastal site; however, year and site had no significant relationship with abundance of sweep-net taxa (Table 2). 2.60 † Secondary guilds always overlapped with a primary guild (e.g., a pollinator was also a predator) so percentages sum to more than 100%.

Overall arthropod biomass
Biomass had a significant positive relationship with willow cover in sweep-net samples but not pitfall samples (Table 3). Neither shrub birch nor ericaceous cover was significantly associated with biomass of either sample type. Slope was the only significant predictor of biomass in pitfall traps, with steeper slopes tending to have less biomass. Biomass did not differ significantly between years or sites for either sample type. Mean biomass varied greatly among taxa and life stages (Appendix S1: Table S2). For pitfall samples, one observation with nearly twice the next highest biomass values was truncated.

Taxon-specific abundance
Most orders were significantly associated with willow cover but not shrub birch or ericaceous cover. Abundance of each order tended to be lower on steep slopes (Appendix S1: Tables S4, S5). Generally, abundance was higher at the coastal site. The effect of year varied, with 2016 having a greater abundance of pitfall taxa but a lower abundance of sweep-net taxa.
We visualized the relationships of these arthropod orders to shrub cover by generating predicted values from the model for each shrub variable while holding all other variables at their mean values (Fig. 3; Appendix S1: Tables S4, S5). We had very few transects with shrub birch cover >25% so model estimates generally had wide prediction intervals for arthropod responses above that level. For Diptera, Hemiptera, Hymenoptera in sweep-net samples, Hymenoptera in pitfall samples, and Lepidoptera (Fig. 3a, c, e, f, and g respectively), abundance increased significantly with willow cover. Hemiptera abundance increased significantly with ericaceous cover, and Hymenoptera abundance in both sample types increased significantly with shrub birch cover. Araneae (Fig. 3b) and Coleoptera (Fig. 3d) were not significantly related to willow cover, but were instead positively (Araneae) or negatively (Coleoptera) associated with ericaceous cover and negatively associated with slope (Appendix S1: Table S4). One observation of Hemiptera fell outside three standard deviations of the mean, despite being a true value, and was thus truncated to the next highest value.

Biomass of trophic groups
Biomass of herbivores was largely hemipteran families, the most dominant being Psyllidae (26.73%) and Cicadellidae (24.58%). Pollinator biomass was largely dipteran families (Culicidae, 28.26%; Empididae, 20.99%; and Chironomidae, 10.95%). Biomass of predators was overwhelmingly Coleoptera (Carabidae, 71.43%; Appendix S1: Table S3). All three major trophic levels were significantly associated with willow cover, but not ericaceous or shrub birch cover ( Fig. 4; Appendix S1: Table S6). The direction of this relationship varied, with high levels of willow cover associated with a lower biomass of predators (Fig. 4a) but a greater biomass of herbivores  (Nakagawa and Schielzeth 2013); R 2 c , conditional pseudo-R 2 (Nakagawa and Schielzeth 2013); SE, standard error. Abundance was analyzed with generalized linear mixed models using a negative binomial distribution. All models included plot as a random effect. Data and model fit are shown in Fig. 2. v www.esajournals.org (Fig. 4c). Pollinator biomass had a significant nonlinear relationship with willow cover, with a predicted peak at about 50% cover and declining at higher levels of willow cover (Fig. 4b). Slope was a significant predictor of predator biomass, with fewer predaceous taxa on steep slopes (Appendix S1: Table S6). Biomass of pollinators was significantly lower in 2016, but we found no Fig. 2. The relationship of willow cover, shrub birch cover, and ericaceous cover to total abundance of arthropods caught in pitfall traps (a-c) and sweep-nets (d-f) in northwestern Alaska, USA. Symbol type indicates site (triangles = coastal, circles = interior) and color indicates year (black = 2015, gray = 2016). Lines (dashed = coastal, solid = interior) show predicted relationship of abundance to each vegetation characteristic from the models described in Table 2. Lines are only shown for significant relationships. v www.esajournals.org 9 May 2021 v Volume 12(5) v Article e03514 significant differences between years for biomass of herbivores or predators. Biomass of trophic groups did not differ significantly between sites.

Community composition
Arthropod community composition was influenced strongly by both willow and ericaceous cover in NMDS analysis. Permutation tests of correlation coefficients indicated that changes in community composition of arthropods caught in pitfall traps were significantly associated with willow cover (R 2 = 0.42, P = 0.001) and ericaceous cover (R 2 = 0.09, P = 0.021), marginally associated with shrub birch cover (R 2 = 0.06, P = 0.086), but not associated with shrub height (Fig. 5). Clustering of plots within each site indicated a high degree of spatial similarity among pitfall arthropod communities (Fig. 5). In sweepnet samples, changes in arthropod community composition were associated with willow cover (R 2 = 0.42, P = 0.001) and ericaceous cover (R 2 = 0.09, P = 0.003) but not shrub birch cover or shrub height (Fig. 6).

DISCUSSION
Our results suggest that expansion of some woody shrubs into tundra-dominated habitats may increase overall abundance of arthropods and the relative abundance of herbivores and pollinators, and therefore, that shrub encroachment has the potential to alter Arctic food webs.
However, shrub taxa differed in their relationships with arthropod abundance, biomass, and community composition. Willow was more strongly associated with changes in community composition and abundance than birch or ericaceous shrubs, supporting our prediction that shrub taxa would vary in their influence on arthropods. Willow, birch, and ericaceous shrubs offer different resources to herbivores and pollinators, and consistent with our prediction, this was reflected in the response of these types of arthropods. Our findings partially corroborate those of previous studies in Arctic Alaska (Rich et al. 2013, in which higher numerical abundance of arthropods was associated with greater shrub cover, a pattern likely driven by greater availability of food as well as structural niches. In contrast to earlier work, we found that arthropod biomass was negatively associated with cover of some shrub taxa. Thus, to obtain accurate predictions of the cascading ecological effects of shrub expansion in a region, one must incorporate not only abiotic variables and dominant plant growth forms but also the identity of dominant plants. Willows (S. alaxensis, S. glauca, S. pulchra, and S. richardsonii) were positively associated with overall abundance, biomass of canopy-dwelling insects, and numerical abundance of most taxa and had the greatest influence of all biotic factors on community composition for ground-dwelling and canopy-dwelling arthropods. It is possible  (Nakagawa and Schielzeth 2013); R 2 c , conditional pseudo-R 2 (Nakagawa and Schielzeth 2013); SE, standard error. Biomass was analyzed with linear mixed models after transforming the response variable (pitfall biomass was cuberoot transformed, sweep-net biomass was log-transformed). All models included plot as a random effect.
†  that there could be plot or regional differences in climate that may drive both the growth of willows and the abundance of arthropod taxa; however, in both the coastal and interior regions we measured a wide range of willow cover and arthropod abundance values that were positively correlated within region. Because shrub taxa seem to differ in their effects, we caution against simply grouping shrubs by growth form during analysis as doing so may obscure important ecological relationships with other organisms (Saccone et al. 2017). Detailed studies of arthropod communities with classifications to genus or even species level would provide further nuance to our understanding of the ongoing shifts in Arctic food web structure and community composition.
We expected pollinators to be associated with ericaceous shrubs and willow, since both are insect-pollinated. In partial support of this idea, pollinators had a positive association with willow cover. Willow and birch shrub thickets contain shade-tolerant forbs (Swanson et al. 1985), which are pollinated by bees (Apidae), small beetles (Cantharidae and Staphylinidae), and some flies (Muscomorpha, Syrphidae;Schmid-Hempel and Speiser 1988, Steinbach and Gottsberger 1994, Lundgren and Olesen 2005, and willow is a source of nectar and pollen (Fox 1992). Contrary to our expectation, however, pollinator biomass was not associated with ericaceous cover, flowering shrubs typically found in tundra habitats. Because both willow thickets and tundra contain resources for pollinators, we expected a difference in pollinator composition between the two habitats but not necessarily a difference in biomass. It should be noted that our arthropod samples were combined across a period encompassing green-up to seed set for all shrub types, and thus flowering of all shrub types occurred during the arthropod sampling period. The greater pollinator biomass in willow thickets could be due to larger body size of pollinators, a shift in composition toward larger taxa, or more pollinators overall. Investigating this pattern further in future research would be valuable. Although we were not able to analyze alder or high levels of shrub birch cover, we postulate that an increase in these wind-pollinated shrubs would result in declines in pollinators unless they also provide habitat for forb species. . Lines are only shown for significant relationships. Envelopes are 95% prediction intervals from the full model. Wide prediction intervals at high levels of shrub cover result from high variance in arthropod response and few observations of shrub cover above 50%.
High levels of willow cover were associated with a greater biomass of herbivores but a lower biomass of predators. Willow is more palatable and less toxic than birch or ericaceous shrubs (MacLean and Jensen 1985, Mulder 1999, Christie et al. 2015, which may explain the positive relationship with herbivorous arthropods. Some ericaceous shrubs (e.g., Vaccinium spp.) are more palatable than others, so one limitation of our study is that grouping ericaceous shrubs may have disguised variations in arthropod responses to different ericaceous shrub taxa. It is possible that these more palatable ericaceous shrubs have a positive relationship with the abundance of arthropod herbivores. Despite the assumption that increased abundance of herbivorous arthropods might create greater prey availability for predaceous taxa, we saw lower predator biomass in areas of high willow cover. We grouped taxa by life stage (e.g., larva, nymph, adult), and different life stages are likely not preferred equally by predators, which could complicate how abundance relates to prey availability. It should be noted that our inferences about predaceous taxa in general are limited, since total predator biomass was overwhelmingly composed of carabid beetles (Coleoptera: Carabidae). Fig. 5. Non-metric multidimensional scaling (NMDS) ordination results for arthropods caught in pitfall traps in northwestern Alaska, USA. Each symbol represents an arthropod community at a single sampling location in a single year. Symbol type indicates site (triangles = coastal, circles = interior), and color indicates plot within site. Axis 1 is primarily associated with willow and shrub birch cover; axis 2 is primarily associated with willow height, ericaceous cover, and ericaceous height. Vectors showing correlation between shrub variables and community composition are overlaid with length corresponding to significance.
Changes in the relative abundance of different trophic groups could have ecosystem-wide consequences. Increased woody vegetation, coupled with higher air temperatures, could result in increased abundance of herbivorous insects and rates of damage to plants. High temperatures have been linked to increased herbivory of dwarf birch in Arctic biomes (Barrio et al. 2017) Increased air temperatures have also been associated with increased rates of herbivory and a greater diversity of foliar damage in fossil records from high-latitude and temperate ecosystems (Currano et al. 2008, Wappler and Denk 2011, Labandeira and Currano 2013. Furthermore, data from the last few decades indicate that warmer air temperatures in forests are associated with increased intensity of insect outbreaks and range extensions of pests (Logan et al. 2003). Shrub expansion, coupled with projected increases in air temperature, may facilitate increased populations of herbivorous insects in the Arctic. We expect the response of arthropod communities to future shrub expansion will be heterogeneous in space because willow, birch, alder, and ericaceous shrubs expand at differing rates across the Arctic (Myers-Smith et al. 2011, Elmendorf et al. 2012. In areas where willow is Fig. 6. Non-metric multidimensional scaling (NMDS) ordination results for arthropods caught via sweep-net in northwestern Alaska, USA. Each symbol represents an arthropod community at a single sampling location in a single year. Symbol type indicates site (triangles = coastal, circles = interior), and color indicates plot within site. Axis 1 is primarily associated with willow cover and ericaceous height; axis 2 is primarily associated with ericaceous cover, shrub birch cover, and shrub birch height. Vectors showing correlation between shrub variables and community composition are overlaid with length corresponding to significance. responsible for most of the increase in woody vegetation, such as western Canada and Arctic Russia (Myers-Smith et al. 2011), we may see more pronounced increases in arthropod abundance and consequently greater benefits to higher-level insectivores in the animal community. Some pollinators may benefit from this increase in willow and associated forbs, but many that are dependent on ericaceous shrubs and other tundra-associated flowering species may be negatively affected. In addition, predaceous taxa were negatively associated with willow cover and may decline with shrub expansion. In areas where wind-pollinated birch and alder shrubs are expanding, such as Arctic Alaska, pollinators may decline. One intriguing pattern in our community composition data is that sweep-net arthropod communities were much more similar across the entire region, whereas pitfall arthropod communities within the coastal area were much more similar to each other than to communities in the interior. Perhaps this reflects differences in dispersal ability, since many taxa caught in pitfall traps are wingless. If so, this suggests that future effects of increased shrub cover on arthropod taxa may have spatially heterogenous effects that depend on local community composition and particularly dispersal ability of abundant taxa. This result should be interpreted with caution, since the greater movement range and dispersal ability of flying insects indicates that they are able to exploit a larger portion of the landscape, and our vegetation sampling (20 m transect) may not have adequately captured the resources available to them across their range.
In addition to spatial heterogeneity, we project that changes in arthropod communities will not be linear over time. We observed a curvilinear relationship of some arthropod groups (Diptera and pollinators overall) with willow cover, with abundance peaking at about 40-60% cover; these results suggest that willow thickets can provide good habitat for these groups unless the shrubs are too dense. We had few observations of shrub cover exceeding 60%, leading to uncertainty in these predictions, but we suspect that dense shrub thickets may be less suitable habitat than areas of low to moderate shrub cover due to reduced solar radiation, reduced plant diversity, and the retention of snow, which can delay green-up, soil thaw, and thus arthropod activity and development (Myers-Smith et al. 2011, Sweet et al. 2014. Shrub thickets provide greater plant biomass and vertical habitat for arthropods and thus tend to host greater abundance, but we have provided evidence that these relationships vary among arthropod and shrub taxa. Increases in Arctic shrub cover may increase food availability for insectivores such as migratory birds (Thompson et al. 2016), which are highly dependent upon the arthropods that flourish in the Arctic during the summer breeding season. Therefore, nonlinear effects are likely to cascade throughout the ecological community.