Parasitoids indicate major climate-induced shifts in arctic communities.

Climatic impacts are especially pronounced in the Arctic, which as a region is warming twice as fast as the rest of the globe. Here, we investigate how mean climatic conditions and rates of climatic change impact parasitoid insect communities in 16 localities across the Arctic. We focus on parasitoids in a widespread habitat, Dryas heathlands, and describe parasitoid community composition in terms of larval host use (i.e., parasitoid use of herbivorous Lepidoptera vs. pollinating Diptera) and functional groups differing in their closeness of host associations (koinobionts vs. idiobionts). Of the latter, we expect idiobionts-as being less fine-tuned to host development-to be generally less tolerant to cold temperatures, since they are confined to attacking hosts pupating and overwintering in relatively exposed locations. To further test our findings, we assess whether similar climatic variables are associated with host abundances in a 22 year time series from Northeast Greenland. We find sites which have experienced a temperature rise in summer while retaining cold winters to be dominated by parasitoids of Lepidoptera, with the reverse being true for the parasitoids of Diptera. The rate of summer temperature rise is further associated with higher levels of herbivory, suggesting higher availability of lepidopteran hosts and changes in ecosystem functioning. We also detect a matching signal over time, as higher summer temperatures, coupled with cold early winter soils, are related to high herbivory by lepidopteran larvae, and to declines in the abundance of dipteran pollinators. Collectively, our results suggest that in parts of the warming Arctic, Dryas is being simultaneously exposed to increased herbivory and reduced pollination. Our findings point to potential drastic and rapid consequences of climate change on multitrophic-level community structure and on ecosystem functioning and highlight the value of collaborative, systematic sampling effort.

INTRODUCTION

Climate change can affect species distributions (Hickling, Roy, Hill, Fox, & Thomas, 2006; Jepsen et al., 2011; Parmesan, 2006; Parmesan & Yohe, 2003), phenology (Høye et al., 2014), and fecundity (Bowden et al., 2015), with knock‐on effects on community composition (Habel et al., 2016; Koltz, Schmidt, & Høye, 2018), on the strength and identity of biotic interactions (Both, van Asch, Bijlsma, van den Burg, & Visser, 2009; Van Nouhuys & Lei, 2004), and ultimately on ecosystem functioning (Ammunét, Kaukoranta, Saikkonen, Repo, & Klemola, 2012; Memmott, Craze, Waser, & Price, 2007; Schmidt et al., 2016; Schmidt, Mosbacher, et al., 2016; Tiusanen, Hebert, Schmidt, & Roslin, 2016). The Arctic region has, on average, been warming about twice as fast as the rest of the globe, making it an important observatory for climate change impacts (IPCC, 2007; Walsh, 2014). The arctic fauna is dominated by arthropods (Høye & Culler, 2018; Schmidt et al., 2017; Wirta et al., 2016), yet temporal data on arthropod diversity and abundance are sparse, limiting current insights into how rapid environmental change is affecting arthropod community structure and functioning (Gillespie et al., 2020). In a rare case where such data are available, a dramatic decline and species turnover have been reported among one of the key pollinating taxa (Diptera: Muscidae; Loboda, Savage, Buddle, Schmidt, & Høye, 2018).

To circumvent the current shortage of arctic arthropod data, we can substitute temporal gradients with geographic ones (Blois, Williams, Fitzpatrick, Jackson, & Ferrier, 2013; Elmendorf et al., 2015; Körner, 2007), comparing regions currently characterized by different climate regimes. The responses of species to gradients over time are likely to mimic their occurrence patterns in space (Romero et al., 2018). This approach does have fundamental limitations in a world where novel communities and environments are formed (Damgaard, 2019). However, in regions of rapid change, recent shifts in climate may already have remodeled insect communities (Rafferty, 2017). This offers an opportunity for a new type of space‐for‐time approach: a comparison among regions experiencing different rates and types of change. (i.e., a space‐for‐change approach). However, few studies to date have exploited such contrasts (but see Prevéy et al., 2017; Scarpitta, Vissault, & Vellend, 2019). In this context, the Arctic offers a rare opportunity, as it simultaneously encompasses strong gradients in climatic conditions and regional variation in the recent rate and mode of climate change, even at a subcontinental scale (Abermann et al., 2017).

In this study, we examine the impacts of both regional mean climate and recent shifts in climate on trophically structured insect communities. As our model system, we use a tritrophic food web, which includes a widespread flowering plant, the mountain avens (genus Dryas in family Rosaceae), its lepidopteran herbivores and dipteran pollinators, and the parasitoids of these herbivores and pollinators. Here, we define parasitoids as predators developing in close association with a single individual of its host species, typically killing it in the process (Hawkins, 2005). We refer to Lepidoptera as herbivores, as the vast majority of lepidopteran species of the Arctic have plant‐feeding larvae, and since we quantify lepidopteran herbivory on our focal plant species (Dryas). We define pollinators as the large array of flower‐visiting Diptera, which provide the main visitors and pollinators of Dryas (Tiusanen et al., 2016, 2019). Our focus on these particular guilds and taxa is motivated by a series of simple considerations. First, insect herbivores are key players in the tundra biome in terms of species richness (Wirta et al., 2016) and biomass (Bar‐On, Phillips, & Milo, 2018). Second, herbivores influence arctic vegetation both through low‐level background damage (Barrio et al., 2017; Rheubottom et al., 2019) and severe episodic and occasional population outbreaks (Lund et al., 2017; Post & Pedersen, 2008), with defoliators such as moth and sawfly larvae causing most of the loss of plant biomass. Third, a large proportion of arctic plants are insect‐pollinated (Kevan, 1972), making plant–pollinator interactions some of the main determinants of arctic insect communities (Tiusanen et al., 2016, 2019). Fourth, both insect herbivore and pollinator populations are at least to some degree regulated by generalist predators, and especially by parasitoids (Letourneau, Jedlicka, Bothwell, & Moreno, 2009). Fifth, parasitoids for their part are expected to be especially sensitive to environmental change due to their high trophic level (Voigt et al., 2003). Based on these considerations, we expect climatic impacts to exert a major effect on arctic parasitoid communities through direct impacts on species and through knock‐on effects mediated by biotic interactions (Schmidt et al., 2017). As groups of parasitoid species show affinities for hosts of certain phylogenetic branches (e.g., Diptera, Lepidoptera, etc.), they are indeed likely to reflect changes in widely different parts of the arctic arthropod food web, effectively serving as sentinels of the arthropod community.

In this paper, we use the structure of parasitoid communities as a surrogate for the structure of the total insect community, and infer the relative abundances of different host taxa from the abundance of different parasitoid taxa for which the hosts are known. This approach is underpinned by a body of literature showing a positive correlation between parasitoid and host abundances across the community (Askew & Shaw, 1986; Godfray, 1994; Hassell, 2000). We explore climatic impacts on two key dimensions of ecological variation among parasitoids: on host‐use taxonomy (i.e., the use of hosts in order Diptera vs. Lepidoptera) and host‐use strategy (Figure 1). To characterize host‐use strategy, we use the simple dichotomy between idiobionts and koinobionts, of which idiobionts kill or paralyze their host at parasitism, whereas koinobionts allow the host to feed and grow during the interaction. Idiobionts often attack non‐growing host stages such as eggs and pupae, whereas koinobiont often attack growing stages (larvae). The koinobiont strategy requires physiological adaptation and therefore restricts the host range (Godfray, 1994). The two groups have been shown to differ also in their overwintering ability in favor of koinobionts. (Hance, van Baaren, Vernon, & Boivin, 2007; for further details, see Section 2.2).

FIGURE 1

Conceptual summary of parasitoid life‐history strategies, host group preferences, and their links to multidecadal mean climate and climate change and its implications. For each of three aspects of parasitoid ecology, that is, parasitoid life‐history strategy, parasitoid host group taxonomy, and associated host abundances, we identify the expected responses to mean climate and recent climate change. We identify the response categories (classes) scored as Attributes contrasted, the biological features of each class as Distinguishing features, and expectations in terms of responses in terms of two types of patterns: changes in the dominance of the respective group with a change in mean conditions (column Mean climate), and changes in the dominance of the respective group with recent trends in a warming Arctic (column Recent change). Finally, we summarize the results obtained in terms of contemporary patterns across the Arctic (Spatial results) and matching patterns in the 22 year time series from Zackenberg, Northeast Greenland (temporal results). Given the dominance of Diptera among arctic pollinators and Lepidoptera among arctic herbivores, we note that changes in host use provide a window to the relative abundance of these key guilds. For clarity, we color code the taxonomically and ecologically separate Lepidoptera‐ and Diptera‐based food web modules in green and blue, respectively, reminding the reader that larval Lepidoptera form the dominant herbivores of Dryas, whereas adult Diptera form the dominant pollinators

Overall, we ask:

How do mean multidecadal climatic conditions affect arctic parasitoid community structure?

We hypothesize that the harsh arctic winter is the main environmental filter dictating insect community composition. If this is the case, then we expect to see regional climatic conditions reflected in the predominance of particular life‐history traits within parasitoid communities. We expect parasitoid species that use dipteran (rather than lepidopteran) larval hosts to be more common toward the High Arctic, mirroring the diversity patterns of the host groups (Böcher, Kristensen, & Pope, 2015). In addition, we expect parasitoid communities of colder sites to be dominated by a koinobiont life‐history strategy, which is associated with increasing specialization and cold‐hardiness. Such a trend in parasitoid community composition with climate would match patterns found across elevation gradients (Maunsell, Kitching, Burwell, & Morris, 2015), as well as previously suggested latitudinal trends among parasitoids (Quicke, 2012; Timms, Schwarzfeld, & Sääksjärvi, 2016) and across mutualistic interactions (Schleuning et al., 2012).

How does recent climate change affect the host use of parasitoid communities?

We expect an increase in dominance of parasitoid species using lepidopteran versus dipteran larvae as hosts in areas where summer temperatures have risen more. Such a prediction is supported by recent findings that warmer and drier conditions in parts of the tundra biome result in dramatic population declines among flies (Loboda et al., 2018), while many lepidopteran species may conversely benefit from warming conditions (Habel et al., 2016; Hunter et al., 2014; Klapwijk, Csóka, Hirka, & Björkman, 2013). We reiterate that arctic pollinators are dominated by adult Diptera (Kevan, 1972; Tiusanen et al., 2016, 2019), whereas larval Lepidoptera are dominant herbivores in our focal Dryas heath habitat. Thus, changes in parasitoid host use will also reflect climatic impacts on the two guilds of pollinators and herbivores, respectively.

How does recent climate change affect life‐history strategies within parasitoid communities?

We hypothesize that recent climatic change in the Arctic has considerably impacted parasitoid community composition, which should thus reflect regional differences in the rate and mode of warming. Specifically, we expect the functional composition of communities in faster warming areas to have shifted more toward that of communities occurring at lower latitudes, that is, an increase in idiobiont strategies. Importantly, this ratio can be examined both with respect to diversity and abundance, where counts of species reflect slower evolutionary and biogeographical processes, and counts of individuals reflect ecological processes. Thus, we expect the effects of recent climate change to be more pronounced in the relative number of idiobiont individuals in the community, but less so in the relative number of idiobiont species, as more generally suggested by the work by Menéndez et al. (2006).

How is the level of herbivory on Dryas associated with mean climate and recent change?

Based on the predictions of shifts in parasitoid community compositions and their life‐history traits (see questions 1–3), we expect a shift toward more idiobiont parasitoid communities to weaken predation pressure on herbivorous Lepidoptera, as idiobionts typically have relatively lower fecundity (Price, 1972). We also expect herbivores to gain direct benefits from warmer temperatures. Both of these processes would translate into increased levels of herbivory (de Sassi et al., 2012).

MATERIALS AND METHODS

Study system

The focal systems within this study were heathlands dominated by species of mountain avens (genus Dryas in the Rosaceae family). Such heathlands are common and widespread throughout the arctic biome (Welker, Molau, Parsons, Robinson, & Wokey, 1997). All three Dryas species (native to the arctic and alpine regions of Europe, Asia, and North America, i.e., D. octopetala, D. integrifolia, and D. drummondi; see Tiusanen et al., 2019) are perennial, cushion‐forming evergreen dwarf shrubs. The flowering phenology and other characteristics of Dryas have been shown to be sensitive to temperature both experimentally (Welker et al., 1997) and through monitoring (Panchen & Gorelick, 2015). Snow and water availability are also key factors modulating the phenological response of Dryas (Bjorkman, Elmendorf, Beamish, Vellend, & Henry, 2015), as is also the nutritional content of the plant (Welker et al., 1997).

Dryas plays a central role in several plant–insect interactions. First, it is a food source to many herbivorous noctuid moth species (Lepidoptera: Noctuidae). In particular, moths of the genus Sympistis are specialized herbivores of Dryas. These moths have a 2 year life cycle, first hibernating as an egg and then as a pupa (Ahola & Silvonen, 2005). Importantly, insect herbivory on Dryas (by Sympistis as well as other species) is concentrated on the flowers, directly affecting the reproductive success of the plants (Figure 1). In this study, we therefore define herbivory as florivory, measured as the proportion of damaged Dryas flowers. Second, a major part of insect taxa (most notably Diptera) within high‐arctic insect communities visit Dryas flowers (Tiusanen et al., 2016, 2019) and subsequently aid its pollination. Dryas has therefore been identified as an interconnecting node at the core of arctic food webs (Schmidt et al., 2017). Because pollination can be directly linked to the reproductive output of plants, any loss of specialist pollinators can have dramatic effects on seed production (Auw, 2007), which applies also to the primarily insect‐pollinated Dryas (Tiusanen et al., 2016).

Both the lepidopteran herbivores and dipteran pollinators of Dryas serve as host species to parasitoid wasps (Hymenoptera) and flies (Diptera: Tachninidae; Várkonyi & Roslin, 2013; Wirta et al., 2015). These parasitoid species may shape the abundance and community composition of herbivores and pollinators, and therefore the performance of Dryas. To clearly distinguish between “parasitic Diptera” (i.e., species in family Tachninidae) and “parasitoids using Diptera as hosts”, we henceforth refer to the former as “dipteran parasitoids” (or “tachinids”) and the latter as “parasitoids of Diptera” or Diptera‐using parasitoids.

Parasitoid biology

Parasitoids are organisms living in close association with a single host individual, which they kill at some stage of their development. The physical association with host species varies, but is broadly categorized into two strategies (Figure 1): idiobionts, which halt development of their hosts, and koinobionts, which allow host development to continue, residing within through successive developmental stages (Askew & Shaw, 1986; Godfray, 1994). These strategies are correlated with a suite of other traits, with ramifications for degree of specialization, potential for top‐down control of host populations, overwintering ecology, and phenology (Quicke, 2015). Compared to koinobionts, idiobionts typically have a wider diet (i.e., attack more diverse host taxa). These parasitoids may not need to track their host as they can rely on other host species. While koinobionts can exploit the overwintering behavior of their mobile hosts, idiobionts are restricted to overwintering either as adults or inside hosts in more exposed habitats (Hance et al., 2007). As a likely consequence, idiobionts are more sensitive to winter conditions (Figure 1), and show more pronounced drop in diversity toward higher latitudes than do koinobionts (Timms et al., 2016).

Beyond the koino‐ versus idiobiont dichotomy, parasitoids vary with respect to their specialization for host phylogeny or lifestyle (Figure 1; Quicke, 2015). Some parasitoid groups are tightly associated with certain host groups like cyclorraphous flies for example, while others prey on the silky structures spun by nearly any arthropod. Parasitoid communities are thus expected to reflect underlying geographical patterns of host taxa. Different groups of insect herbivores show distinct diversity patterns, with Lepidoptera (moths and butterflies) dominating at low latitudes (Kerr, Vincent, & Currie, 1998) and Symphyta (sawflies) reaching their highest diversity toward higher latitudes (Kouki, Niemelä, & Viitasaari, 1994). The same applies for other functional groups, such as pollinators. In the High Arctic, pollinator communities predominantly consist of Diptera, while at lower latitudes, hymenopteran pollinators (bees) are regarded as the most important (Böcher et al., 2015). A conceptual summary of parasitoid strategies and host associations and their spatial and temporal implications is presented in Figure 1.

Empirical data

To resolve how climate shapes plant–insect–parasitoid interactions across the Arctic, we set up a distributed, standardized sampling design (Figure 2). We used molecular species identification to characterize parasitoid communities, and globally available remote‐sensed climate data to examine how these communities are structured relative to long‐term abiotic conditions and to recent climate change. To assess the links between food web structure and ecosystem functioning, we quantified herbivory on a key plant resource (larval feeding damage on flowers in the genus Dryas). We then used a local time series of plants, herbivory, and insects to confirm that the predictions derived from large‐scale patterns of parasitoid community composition and herbivory are also observable through time.

FIGURE 2

The structure of the dataset and the links between data sources. The box on the left summarizes data collected across the Arctic on parasitoid community composition and level of herbivory. Parasitoid communities were characterized by host use and parasitoid life‐history strategy (as nested within host use). These spatial data were collected at each of 19 field sites, identified by pink markers on the central map. For each of these sites, we also extracted two types of climate data: variables describing mean temperature and precipitation over the time period 1970–2000 (illustrated in upper hemispheres) and variables describing the rate of the recent temperature change during 2000–2017 (illustrated in lower hemisphere). The box on the right summarizes data used to analyze temporal patterns of host availability at one of the study locations (Zackenberg, Northeast Greenland). The data encompass local climatic data since 1996–2017, counts of muscid flies in insect traps, and annual peak fractions of damaged Dryas flowers on permanent monitoring plots. For clarity and consistency with Figure 1, we show parasitoids in black, pollinators in blue and herbivores in green. Numbers identify sampling localities: 1. Zackenberg, 2. Churchill, 3. Igloolik, 4. Bylot Island, 5. Qeqertarsuaq/Disko Island (low and high altitudes), 6. Kangerlussuaq (low and high altitudes), 7. Kangerluarsunnguaq/Kobbefjord, 8. Hochstetter Forland, 9. Snæfellsnes, 10. Ny‐Ålesund, 11. Svare/Vågå, 12. Finse, 13. Kevo, 14. Finnmark (two different mountains), 15. Monchegorsk, 16. Yamal. For detailed site‐specific information, see Table S1. We note that data from the Russian‐Canadian Arctic are very sparse, reflecting logistic challenges during the focal study period (summer of 2016). For consistency with Figure 1, we color code the taxonomically and ecologically separate Lepidoptera‐ and Diptera‐based food web modules in green and blue, respectively, reminding the reader that larval Lepidoptera form the dominant herbivores of Dryas, whereas adult Diptera form the dominant pollinators

Distributed study design

To measure herbivory and sample parasitoid communities across the Arctic, we joined together through the arctic research networks INTERACT (https://eu‐interact.org/) and NEAT (http://neat.au.dk/). The former is a large‐scale consortium specifically aimed at tackling large questions in arctic research by drawing on its 70+ member stations, and the latter is a recent collaboration of arctic arthropod ecologists. Before the summer of 2016, we sent comprehensive sampling sets to teams in 21 field locations (Figures 2 and 3; see http://www.helsinki.fi/foodwebs/parasitoids/interactsamplinghires.pdf). This sampling was either cancelled or relocated in three of the intended locations due to unexpectedly early spring phenology in Alaska, western Canada, and southernmost Greenland. Thus, the dataset comprises of data from 19 sampling localities, all of which are used in the analyses of community structure and 16 in the analysis of herbivore damage (see Supporting Information [Link], [Link] for details).

FIGURE 3

Relationship between host use (y axes), the rate of temperature change in the winter period (x axes) and in the summer period (in panel c), with the three curves corresponding to the models estimates for low, mean, and high values occurring in the data as indicated in the right‐hand side box, and the colored areas around them showing 95% confidence intervals. The colors of the data points show the local rate of temperature change for the summer period, adhering to the color scheme of the left‐hand legend. Panels (a) and (b) show the effects of these variables on the fraction of parasitoid species and individuals, respectively, which mainly use lepidopteran hosts. Panels (c) and (d) visualize the same trends but for parasitoids of Diptera. The size of the data points is proportional to the number of species or individuals, respectively whereas the colors of data points represent local rate of temperature change for the summer period

Participants were instructed to begin sampling when Dryas started flowering widely in the landscape. Each participant set up three to four sampling plots of 10 white sticky traps (4.5 × 5 cm each), covered by individual wire cages to exclude bird predation on insect catches. Depending on the abundance of Dryas flowers, three to five of these cages acted as landmarks for circular plant monitoring sub‐plots with a radius of 27 cm (52 cm in Zackenberg). Within these subplots, participants counted and scored Dryas inflorescences in four different developmental stages: dark buds, buds with visible petals, open flowers, and senescent flowers. Participants also counted the number of flowers damaged by insect herbivores. During three visits, typically interspersed by 6 days, the sticky traps were set up, changed, and collected, and plant data were collected. Thus, in total, 351 plant plots were surveyed three times each, with a sum of 1,500 sticky traps deployed.

Climatic variables of geographic scale

To characterize multidecadal mean environmental conditions of each sampling site (henceforth “mean climate”), we used climatic data extracted from the online database WorldClim (Fick & Hijmans, 2017) at 0.05° (1.9 × 5.6 km at 70°N) resolution, based on the coordinates of the field collections, including a 0.05° buffer zone. These data were used to characterize the mean climatic conditions (1970–2000; Figure 2), which naturally cover a wide range of values as our study sites span from subarctic to high‐arctic conditions. For a validation of these data, see Figure S3.

To describe recent changes in climate, we calculated linear temperature trends for the past 18 years. For this, we used the monthly average daytime land‐surface temperature supplied by the MODIS satellite platform at a spatial resolution of 0.05° (https://doi.org/10.5067/MODIS/MOD11B3.006; Wan, 2014; Wan, Zhang, Zhang, & Li, 2004), from which we extracted the values including a 0.05° buffer zone around sampling sites. These data describe the changes in land‐surface temperature between 2000 and 2017 (Figure 2). To find clues on possible mechanistic underpinnings of observed patterns, we focus on the summer period (June–August), during which insect reproduction takes place. To account for climatic variation during insect overwintering, we also use the rate of surface temperature change during the winter period (September–May).

Due to the arctic amplification, temperature rise is higher at higher latitudes. In our dataset, this is evident especially in the summer period temperatures, which have risen more at localities with colder mean summer temperatures. Winter temperature changes show a more idiosyncratic pattern with a hot spot region around the Barents Sea. The rates of temperature change during different parts of the year are often correlated at the site level. Surprisingly, the recent temperature change of the summer and the autumn periods was negatively correlated with each other (r = −.74). This correlation is considerably reduced when considering the whole winter period (for a summary of all explanatory variables used during modeling, see Table S3).

To validate that the remote‐sensed metrics used in the analyses are actually reflective of local conditions, we compared our metrics to ground‐level measurements where available (Figures S3 and S4), and compared multiple metrics of recent climate change to each other (Figure S4). Many of these metrics proved highly correlated, providing evidence that they provide consistent and biologically relevant measures of local conditions.

Finally, to control for effects of weather during sampling, we recorded temperature in situ during the trapping period using EL‐USB‐2 temperature loggers (Lascar Electronics) exposed close to the soil surface (see Supporting Information [Link], [Link]).

Sample management and molecular workflow

Samples of insects on sticky traps were stored at −20°C until used for DNA extraction. We then used orange oil (Romax Glue Solvent; Barrettine) to dissolve the glue, moving parasitoid wasps and flies individually onto 96 deep well lysis plates for DNA extraction (NucleoSpin® 96 Tissue kit, REF 740741.4; Macherey‐Nagel; Data1; for a summary of individual datasets, see Supporting Information [Link], [Link], section “Sequence data processing”). The same traps were subsequently rechecked for parasitoids, and additional samples processed using a salt extraction protocol following Kaunisto, Roslin, Sääksjärvi, and Vesterinen (2017; Data 2; for details, see Supporting Information [Link], [Link]). These two DNA extraction methods yielded comparable results (success rates in final data 91.1% and 94.9%, respectively).

To identify the parasitoids, we used primer pair BF and HCO2198, as targeting a variable region of the mitochondrial COI (cytochrome oxidase 1) gene (Table S2). All samples were sequenced using Illumina technology, with the exact workflow and bioinformatics pipeline identified in the [Link], [Link] (Figure S2). Sequences were assigned to operational taxonomic units (OTUs; henceforth “species”) with at least family‐level taxonomic affinities, and the site‐by‐taxon data were used for downstream analyses of taxonomic diversity, biotic niches, and functional composition.

Scoring of parasitoid life‐history traits

Since our sampling sites shared relatively few species, species‐level analysis is less informative for detecting climatic impacts on community‐level patterns. Instead, we describe the functional community composition as (a) the prevalence of parasitoid taxa with different parasitism strategies (idiobionts vs. koinobionts) and (b) the prevalence of parasitoid taxa using different main host groups (Lepidoptera vs. Diptera; Figures 1 and 2). To classify parasitoids by their host‐use and parasitism strategy, we gathered information on these traits from the literature (Böcher et al., 2015; Quicke, 2015; Stireman, O'Hara, & Wood, 2006; Timms, Bennett, Buddle, & Wheeler, 2013; Várkonyi & Roslin, 2013; Yu, Van Achterberg, & Horstmann, 2005). Parasitoids of taxa other than Lepidoptera and Diptera, such as Araneida, Coleoptera, Hemiptera (aphids in particular), Symphyta, or other parasitoids, were not analyzed separately (i.e., as was done for parasitoids of Diptera and Lepidoptera) since their numbers were too low—accounting for 0.7%–7.1% (mean 3.3%) of species and 0.03%–4.3% (mean 1.9%) of individuals per host order. A bigger group left outside of the two focal diet categories were those parasitoids with very wide diets or those that use various hosts within the taxonomic resolution of our identification, accounting for 31.5% of species and 20.1% of individuals. To make our response variables representative of the full parasitoid community, parasitoids of taxa other than Lepidoptera and Diptera were still included in the denominator of our response variables (i.e., within the totals of parasitoid species and individuals, respectively). The sources and criteria used in the trait classifications are further specified in the [Link], [Link] (Appendix S4).

Temporal data

The pan‐arctic data described above provide a single view of current patterns, all derived within a single year. To supplement this snapshot with data on temporally resolved changes within a particular site (Figure 2), we used herbivory and arthropod abundance data collected from 1996 to 2017 at Zackenberg (74°28ʹN, 20°34ʹW) by the BioBasis monitoring program (Schmidt, Hansen, et al., 2016; Schmidt, Mosbacher, et al., 2016). These data were provided by ZERO (Zackenberg Ecological Research Operations) and Asiaq—Greenland Survey and are available in the GEM database (http://g‐e‐m.dk/) and summarized in annual reports of the Zackenberg Research Station (http://g‐e‐m.dk/gem‐localities/zackenberg/publications/annual‐reports/). For each year, they encompass standardized observations of herbivore damage, insect abundance, and associated environmental variables, including site‐specific snow cover at six regularly monitored plots originally selected to represent different snow conditions. Weekly observations from each plot provide a detailed description of herbivory on Dryas, including counts of buds, fresh flowers, and senescent flowers, and feeding marks by herbivores. For levels of herbivory, we used the annual peak herbivory rates (across the three sampling events) reported annually.

To quantify changes in the dominant host taxa available to parasitoids in the Zackenberg insect community, and of pollinators available to plants, we extracted data on the abundance of muscid flies (Diptera: Muscidae), measured as the number of individuals caught in two window traps and 20 pitfall traps each summer before 26th of August (data available at GEM database; Figure 2).

As potential determinants of the level of flower damage and pollinator abundance, we used the date of snow melt available at the plot level, and yearly values of summer and autumn soil minimum temperatures at 10 cm depth (measured at the nearby climate station). Here, we note that the monitoring plots were originally chosen to represent an environmental gradient from snow‐accumulating depressions to windswept areas (Schmidt, Hansen, et al., 2016; Schmidt, Mosbacher, et al., 2016), and that they thus vary substantially in the relative timing of snow melt. To characterize local temperature conditions, we use soil temperatures, since they capture the summer temperatures as experienced by insects. By integrating the combined effect of ambient temperatures, solar warming, and water content, these temperature data are akin to the surface infrared reflectance captured by satellites in our large‐scale data (see climatic variables above). For the autumn period, soil temperatures record overwintering temperatures experienced by the insects due to combined effects of ambient temperatures and the presence or absence of snow cover. To fully account for the temperatures experienced by the insects, we also included the focal summer ambient air temperature (at 2 m height), which is the most important factor in explaining the activity of flying insects within the season (Høye & Forchhammer, 2008).

Statistical analyses

To quantify the effects of mean climate and recent climate change on functional descriptors of parasitoid community composition, we used generalized linear modeling, with the fraction of parasitoids adhering to a given main host (Diptera or Lepidoptera) or host‐use strategy (prevalence of the idiobiont strategy) as our response variables (see Figure 1). As we expected differential responses at the level of species and individuals (see Section 1), we modeled the fraction of individuals in the community and the fraction of species in the community as separate responses. All in all, we modeled eight response variables (Table S4). The first four of these reflect parasitoid host use: the fraction of species (model 1A) and individuals (model 1B) of parasitoids associated with Lepidoptera, and the fraction of species (model 2A) and individuals (model 2B) of parasitoids associated with dipteran hosts. The last four models focus on the parasitism strategy, reflecting the fraction of parasitoid species (models 3A and 4A) and individuals (models 3B and 4B) adhering to an idiobiont strategy for parasitoids of Lepidoptera and Diptera, respectively.

In modeling each of these response variables, our overarching objective was to evaluate the evidence for imprints of both baseline climatic conditions and recent change. Since climatic impacts may relate to multiple different climatic descriptors as calculated for several parts of the year, there is a nontrivial risk of overfitting. To this aim, we placed special emphasis on the logic and sequence of model building, as further explained in Appendix S9. In brief, we first used univariate analysis to assess individual explanatory variables describing multidecadal mean climatic conditions independently from each other. We then tested if adding a second variable describing mean climatic conditions improved the fit. Finally, we tested for added effects of recent climate change. These variable families were then brought into a joint model in a hierarchical manner, starting from the longer term impact (averages of winter temperature, summer temperature, winter precipitation, and summer precipitation) and progressing to recent change (rate temperature change in winter and summer). Variables were retained for the final model based on the reduction in QAICc observed (Lebreton, Burnham, Clobert, & Anderson, 1992). QAICc is a quasi‐likelihood counterpart to the corrected Akaike information criterion (AICc), and is better suited for modelling overdispersed count or binary data. QAICc values were calculated with R package MuMIn (Barton, 2016), calculating a global dispersion parameter from a model containing all of the variables included in the models being compared at a time. By this overall approach, we specifically answer the following chain of questions: First, do long‐term conditions affect the focal community descriptor (response)? Second, with these impacts accounted for, do metrics of recent change add additional explanatory power? The model selection process is summarized in Table S5.

To model the level of herbivory in sampling plots across arctic localities, we fitted a binomial mixed‐effects model (model 5) with the maximum fraction of herbivore‐damaged flowers recorded for each sampling subplot as the dependent variable. We used the same selection of fixed climatic effects as in the models of functional parasitoid community composition (see Table S3) and constructed the model sequentially starting from mean climate variables and subsequently testing if variables describing climatic change improve the overall fit. Since it is difficult to distinguish damage from senescent flowers, and large altitude differences between plots is bound to affect unrecorded local conditions, we included the mean phenological phase when plant surveys were done and the relative altitude within a locality as additional fixed effects. Furthermore, since the abundance of Dryas might affect either the presence of specialist herbivores or saturate damage in dense flower stands, we also included the logarithm of flowers recorded for each subplot as a fixed effect. We introduced these three methodological variables as the null model prior to the sequential model construction. To account for random variance in the intercept between sites and plots, we introduced plots within localities as a random effect. To correct for overdispersion, we also included a random intercept effect at the observation level (random residual). The model selection process is summarized in Table S6. To facilitate the interpretation of the results, we report the effect sizes as odds ratios (OR), which are calculated by exponentiation of the linear predictor. In other words, the OR identifies the change in the odds of the modeled outcome for each unit increase of the explanatory variable (1 SD in the case of standardized variables; Rita & Komonen, 2008).

To evaluate whether the inferred drivers of flower herbivory patterns in space (model 5) also generate similar patterns in time, we used the 22 year time series of Dryas damage from Zackenberg, Greenland. Here, the fraction of Dryas flowers damaged per year in each of the six monitoring plots was modeled by logistic regression (model 6), using matching climatic descriptors as above (see Section 2.3.5 above): mean air and soil temperatures for summer months (June–August) and mean soil temperature during the previous autumn (September–November). To account for the 2 year development of the most abundant herbivore species (Sympistis zetterstedtii (Staudinger, 1857), Lepidoptera: Noctuidae), we included time‐lagged versions of the summer and autumn soil temperatures, shifted either by a year or two. The autocorrelative effect of the response variable was also modeled by including time‐lagged counts of damaged flowers. Finally, we included the relative snow melt date of the monitoring plot in a given year (see Section 2.3.5 above). To account for repeated measures from the same plot, we introduced a random intercept at the level of the monitoring plot, and a random residual effect at the level of observations, to account for overdispersion.

Finally, to evaluate how the same climatic variables as above affected the abundance of key pollinators (muscid flies) over time, we fitted a generalized mixed‐effects model to the number of flies caught per season per trapping station (model 7). As explanatory fixed factors we used the mean air and soil temperatures for the summer months (June–August), the mean soil temperature during the previous autumn (or early high‐arctic winter; September–November), and the 1 and 2 year time‐lagged versions of summer and autumn soil temperatures (used to catch carry‐over effects of soil moisture and population trends). In addition to temperature variables, we included the relative snow melt date on the trapping stations, the trap type (window traps at one trapping station compared to pitfall traps at the other four), and the total number of trapping days, which varied slightly between years and trapping stations. To account for repeated measures from the same site, we introduced a random intercept at the level of the trapping station. To scale the residual variance to match a Poisson error distribution, we included an observation‐level random intercept effect. To facilitate the interpretation of the results, we also report the effect sizes by exponentiation of the linear predictors, which gives the multiplicative change with every unit (here: with each standard deviation) increase in the explanatory variable.

All mixed‐effects models (summarized in Table S4) were fitted with R (R Core Team, 2019) package lme4 (Bates, Mächler, Bolker, & Walker, 2014). We note that the modeling approach followed a clear‐cut logical structure: for explanatory variables measured at the level of the sampling locality only (models 1–4), we modeled data at the site level only. Models 5–7 concern fractions of units observed sharing a particular fate. Model 5 includes plot‐level flower numbers, possible local altitudinal gradients, and senescence data, which is why we explicitly modeled effects at both hierarchical levels (site and plot). Models 6 and 7 focus on data of an entirely different structure, as they were adopted from the Zackenberg time series data, with multiple hierarchical levels.

RESULTS

Altogether, we collected 6,009 parasitoid specimens from the 19 sampling locations across the Arctic (Figure S1a; Table S1; Supporting Information [Link], [Link]). Locations differed greatly in the number of parasitoid individuals caught during the trapping period (Figure S1b; Table S1; Supporting Information [Link], [Link]), a pattern partially explained by local weather conditions during the specific sampling days (Tables S1, S7, and S8). The success rates of parasitoid identification from molecular data were high (Table S1), with 93% of a total 4,699 samples yielding an identifiable DNA barcode sequence. In this material, we detected 460 parasitoid OTUs, of which 80% (90% of successfully barcoded individuals) were attributable to a named genus. The proportion of taxa identified to levels above species did not vary systematically among sites.

Climatic impacts on parasitoid community composition

We detected strong impacts of both mean climatic conditions and recent climatic change on the composition of parasitoid communities across the Arctic (for a summary, see Figure 1).

The colder the mean summer temperature at a site was, the larger was the proportion of Lepidoptera‐using parasitoid species in its parasitoid species pools. On top of this, there was a stronger effect of climate change than of mean climate: parasitoid communities at sites experiencing recent warming in winter period demonstrated a low relative abundance of parasitoids of Lepidoptera. This pattern was evident at the level of parasitoid species (OTUs; Figure 3a; Table 1: M1a) and suggested for individuals (Figure 3b; Table 1: M1b). Conversely, parasitoids of Diptera showed the opposite trend, with a larger fraction of Diptera‐using parasitoid species in areas where winter temperatures have risen the most, while rate of temperature change in the summer period had a negative effect on the proportion of Diptera‐using parasitoid species (Figure 3c; Table 1: M2a). At the level of individuals, the rate of temperature change in the winter had the single largest explanatory effect on Diptera use (Figure 3d; Table 1: M2b). Both effect size and uncertainty were higher at the level of individuals than species.

TABLE 1

Summary of models of functional parasitoid community composition (with explanatory variables defined in Table S3 and models identified in Table S4). Rows show the full set of variables considered in model selection, whereas cell values identify estimates for terms retained on the basis of their QAICc values, with parameter estimates from the resulting, final model (see Section 2 for details). For this table, covariates have been standardized to a mean of 0 and an SD of 1. The statistical significance of intercept and slope estimates are indicated by asterisks, † p < .1 *p < .05, **p < .01, ***p < .001, with significant values (p < .05) highlighted in bold face

Response variable as a fraction of:	Intercept		Explanatory variable
			Winter temperature		Summer temperature		Winter precipitation		Summer precipitation		Winter temp. change		Summer temp. change
	Estimate	SE	Estimate	SE	Estimate	SE	Estimate	SE	Estimate	SE	Estimate	SE	Estimate	SE
Parasitoid species of Lepidoptera (M1a)	−1.19***	0.17			−0.46*	0.17					−0.52*	0.23	0.39†	0.19
Parasitoids of Lepidoptera (M1b)	−1.66**	0.45									−1.18	0.75
Parasitoid species of Diptera (M2a)	−1.06***	0.08									0.18*	0.07	−0.33**	0.09
Parasitoids of Diptera (M2b)	−0.22	0.29									0.63*	0.25
Idiobiont species of parasitoids of Lepidoptera (M3a)	−0.80***	0.16	0.68**	0.20									0.41*	0.19
Idiobiont individuals of parasitoids of Lepidoptera (M3b)	−0.86**	0.24	1.67***	0.32			−0.50†	0.28					0.58*	0.24
Idiobiont species of parasitoids of Diptera (M4a)
Idiobiont individuals of parasitoids of Diptera (M4b)	−1.39***	0.21			−1.27***	0.25			0.62*	0.26	−0.38**	0.11

In terms of parasitism strategies, the fraction of idiobiont species varied as hypothesized (question 3) among the parasitoids of Lepidoptera. On average, communities at warmer climates and at lower latitudes were characterized by a higher fraction of idiobionts. Yet, on top of this trend, areas that had experienced stronger increases in summer temperatures held more idiobiont species than expected based on their mean climate (Figure 4a; Table 1: M3a). At the level of individuals, the fraction of idiobionts was still restricted by low minimum temperatures, but increased more steeply with faster summer temperature changes than did the number of species (Figure 4b; Table 1: M3a). This suggests that a recent change in species abundances rather than a long‐standing status quo is responsible for the observed pattern. As for the species proportions among parasitoids of Diptera, we detected no significant impact of climate on the fraction of idiobionts (Table 1: M4a). In the total pool of Diptera‐using parasitoid individuals, the dominance of idiobionts quickly diminished toward warmer mean temperatures and areas where winter temperatures showed higher rates of increase.

FIGURE 4

Relationship between the functional community composition of the parasitoids of Lepidoptera as the fraction of idiobionts (y axes), the average of mean winter temperatures (x axes), and the rate of change in summer temperatures (with the three curves corresponding to the models estimates for low, mean, and high values occurring in the data, and the colored areas around them showing 95% confidence intervals). Panel (a) shows the model‐fitted effects of these variables on the fraction of idiobionts out of all species of primary parasitoids of Lepidoptera and panel (b) shows the same relationship, but for the fraction of idiobionts out of all individuals of primary parasitoids of Lepidoptera. The size of each data point is proportional to the number of (a) species or (b) individuals, respectively. The colors of data points represent the rate of summer temperature change at the respective locality (see legend on the right)

Herbivory levels across the Arctic

Across the Arctic, the fraction of Dryas flowers damaged by lepidopteran larvae varied substantially but was more closely associated with patterns of recent climate change than with the mean regional climate (for a summary, see Figure 1). A part of the variation observed in herbivory was attributable to the rate of summer temperature change (Table 2; Figure 5), while accounting for the effects of the altitude of the subplots and/or their abundance of flowers. Uncertainty with respect to the exact effect of summer temperature change was considerable (OR of 1 SD step [OR] 9.21, confidence limits 2–41).

TABLE 2

Factors affecting the fraction of Dryas flowers damaged by herbivores across arctic sites (Model 5; see Tables S3 and S4). Shown are coefficient estimates, standard errors, and 95% credible intervals for fixed effects. Rows show the full set of variables considered in model selection, whereas cell values identify estimates for terms retained on the basis of their QAICc values, with parameter estimates from the resulting, final model (see Section 2 for details). For this table, the values of explanatory variables have been standardized to a mean of 0 and an SD of 1. Variables for which no values are shown were not retained during model selection. The statistical significance of intercept and slope estimates is given as p value, with significant values (p < .05) highlighted in bold face

Response: Fraction of Dryas flowers eaten	Estimate	SE	95% CI		p value	OR	OR 95% CI
Covariate			Lower	Upper			Lower	Upper
Intercept	−4.59	0.57	−5.71	−3.48	<.0001
Winter temperature	0.15	0.67	−1.17	1.47	.828	1.16	0.31	4.34
Summer temperature
Winter precipitation
Summer precipitation
Winter temperature change
Summer temperature change	2.22	0.77	0.71	3.73	.004	9.21	2.03	41.74
Mean percentage of senescent flowers	−0.37	0.25	−0.86	0.13	.151	0.69	0.42	1.14
Altitude difference within from locality mean	−0.37	0.25	−0.86	0.13	.269	0.69	0.42	1.14
Log(flowers in the plot)	−0.82	0.26	−1.32	−0.31	.001	0.44	0.27	0.73

FIGURE 5

The relationship between the fraction of flowers damaged by lepidopteran herbivores in Dryas plots and the rate of summer temperature change across arctic localities. The size of the marker illustrates the number of Dryas flowers in the survey plot. Color shades illustrate overlapping data points.

Temporal patterns of lepidopteran and dipteran host availability

Large‐scale impacts of climatic variation in space were matched by impacts of year‐to‐year variation in the Zackenberg time series (for a summary of patterns detected, see Figure 1). Here, Dryas damage by lepidopteran larvae increased significantly during summers with warm air temperatures (OR 1.82) and warm (dry) soils (OR 1.44), which were preceded by warm summers (OR 2.72) and cold (snow‐free) autumns affecting the previous generation of moths (OR 0.28; Table 3; Figure 6). Furthermore, the level of flower damage was higher in early‐melting, exposed plots (OR 0.51). The previous years' abundance of moth larvae had no detectable effect on the focal year's damage level.

TABLE 3

Factors affecting the fraction of Dryas flowers damaged by herbivores in the Zackenberg time series (Model 6; Tables S3 and S4). Shown are coefficient estimates, standard errors of those estimates, 95% confidence intervals, and p values for fixed effects. To facilitate interpretation, estimates at the logit scale are also converted to odds ratios (OR) and associated confidence intervals. For this table, variable values have been standardized to a mean of 0 and an SD of 1. Variables for which no values are shown were dropped during model reduction. The statistical significance of intercept and slope estimates is given as p value, with significant values (p < .05) highlighted in bold face

Covariate	Estimate	SE	95% CI		p value	OR	OR 95% CI
			Lower	Upper			Lower	Upper
Intercept	−4.14	0.23	−4.59	−3.69	<.0001
Summer air temperature	0.60	0.17	0.27	0.93	.0003	1.82	1.31	2.54
Summer soil temperature	0.37	0.16	0.04	0.69	.03	1.44	1.04	1.99
Summer soil temperature t‐1	1.00	0.17	0.67	1.33	<.0001	2.72	1.95	3.79
Summer soil temperature t‐2
Previous autumn soil temperature
Previous autumn soil temperature t‐1
Previous autumn soil temperature t‐2	−1.27	0.16	−1.59	−0.95	<.0001	0.28	0.20	0.39
Relative timing of snowmelt	−0.67	0.22	−1.09	−0.25	.002	0.51	0.34	0.78
Flower damage t‐1
Flower damage t‐2

FIGURE 6

Temporal patterns in herbivory and pollinator abundances as observed at Zackenberg, Northeast Greenland. Panel a) shows chronological patterns in the level Dryas damage by Sympistis larvae and the abundance of muscid flies caught at Zackenberg. The solid lines show the actual mean peak percentage of damage recorded and the mean number of muscid flies caught in a trapping station during a summer season. The shaded areas show the confidence intervals of fitted values from models 6 and 7, respectively. The Dryas damage by lepidopteran larvae is shown separately for early plots (light green) and late plots (dark green). For comparison, surfaces in panels (b) and (c) illustrate the effects of two explanatory climatic variables shared between the two models: the air temperature during the focal summer and the soil temperature of the summer 2 years earlier, for Dryas damage and muscid fly abundance, respectively. Note that in panel (a), there is a gap in the line for muscid flies at year 2010. In this year, all arthropod samples were unfortunately and mysteriously lost in transit between Zackenberg and Aarhus, before being sorted, counted, or databased

The same environmental variables that were associated with high lepidopteran herbivore abundance in the Zackenberg time series were also related to decreases in the abundance of muscid flies (Figure 6c). Warmer air temperatures in the current summer were associated with low abundance of muscid flies, with a 13.4% reduction per centigrade. Also opposite to the patterns found for herbivory, muscid flies benefitted from consecutive years with warm autumn soils (i.e., soils insulated by snow: 10.2% and 10.5% per °C increase for the autumn preceding the previous year and the year before that, respectively; Table 4; Figure 6). Against expectation, a higher number of trapping days were associated with fewer individuals caught per season, with a 4% decrease for every additional week of trapping (28 trap days). The trends observed in muscid fly abundances were the same in both window and pitfall traps, suggesting that abundance rather than activity (e.g., the frequency with which flies walk on the ground) was effectively measured.

TABLE 4

Factors affecting the number of muscid flies caught in yellow pitfalls across time at Zackenberg (Model 7; Tables S3 and S4). Shown are coefficient estimates, standard errors of those estimates, 95% confidence intervals, and p values for fixed effects. To facilitate interpretation, estimates at the log‐scale are exponentiated and associated confidence intervals. For this table, variable values have been standardized to a mean of 0 and an SD of 1. Variables for which no values are shown were dropped during model reduction. The statistical significance of intercept and slope estimates is given as p value, with significant values (p < .05) highlighted in bold face

Covariate	Estimate	SE	95% CI		p value	eβ	OR 95% CI
			Lower	Upper			Lower	Upper
Intercept	6.38	0.05	6.28	6.48	<.0001
Summer air temperature	−0.13	0.05	−0.22	−0.04	.006	0.88	0.80	0.96
Summer soil temperature
Summer soil temperature t‐1	−0.08	0.06	−0.19	0.03	.13	0.92	0.82	1.03
Summer soil temperature t‐2
Previous autumn soil temperature
Previous autumn soil temperature t‐1	0.19	0.06	0.08	0.30	.0009	1.21	1.08	1.36
Previous autumn soil temperature t‐2	0.18	0.05	0.07	0.28	.0007	1.19	1.08	1.32
Relative timing of snowmelt
Number of trap‐days	−0.09	0.05	−0.18	0.00	.043	0.91	0.83	1.00
Trap type

DISCUSSION

In this study, we found a distinct imprint of climatic conditions on the parasitoid communities of the Arctic. A priori, we had hypothesized that the harsh arctic winter would be the main environmental filter dictating insect community composition (Section 1; Figure 1). In terms of host use, we found our hypothesis to be too simplistic. In terms of lifecycle strategies, warmer localities across the Arctic were characterized by parasitoid communities with a higher prevalence of the idiobiont parasitoid strategy, thus supporting the pattern found by Timms et al. (2016)—at least for parasitoids of Lepidoptera. Our second hypothesis, that recent climate change would already have affected parasitoid community composition proportionately to the magnitude and mode of the change, was indeed supported: We found that localities which have experienced a faster rise in summer temperatures are currently dominated by parasitoids of Lepidoptera (question 2), exhibit a higher prevalence of idiobionts (question 3), and experience higher levels of herbivory (question 4; for a summary see Figure 1).We also detected a similar signal over time: in a 22 year time series of herbivore damage and pollinator abundance in North‐East Greenland, we found that higher summer temperatures, coupled with cold autumn soils, are linked with increased herbivory by lepidopteran larvae, and with declines in the abundance of dipteran pollinators. This suggests that rapid environmental change in the Arctic is increasing herbivory and reducing pollination of Dryas, a widespread plant species utilized by many arthropods. However, the effect of the total climate change on insects depends on the distribution of warming effects between the summer and winter periods.

Rapid summer warming changes prevailing host use

Our results demonstrate the impact of climate change on parasitoid communities through host associations, as patterns of recent changes override the much larger differences in the regional mean climate (question 1 vs. question 2). At the level of parasitoid species numbers (OTU diversity), winters remaining cold coupled with warming summers were associated with a higher proportion of lepidopteran host use (Figure 3a), while the opposite was found for parasitoid use of Diptera (Figure 3c). At the level of parasitoid individuals (abundance), the effects proved more pronounced, with steeper slopes predicting communities consisting mainly of parasitoids of either Lepidoptera or Diptera at the extremes of our data range (Figure 3b,d). Our estimates (slopes) of such effects came with substantial uncertainty (i.e., wide confidence limits). One reason for this is likely the noise caused by the fact that we used a single proportion to characterize the local dominance of a given host use (individuals representing the focal parasitoid group out of all parasitoid individuals). In reality, this summary fraction is composed of individuals of several different taxa, which may differ substantially in species‐specific abundances. Our focal parasitoid communities also include parasitoids of, for example, aphids or fungus gnats. Since such parasitoids can track the abundances of their hosts even within a season (Mukai & Kitajima, 2019; Nakata, 1995), they may occasionally appear in great numbers (T. Kankaanpää, personal observation), thereby adding noise to local estimates.

How host‐use composition ties into local community structure is illustrated by the time series collected at Zackenberg, Greenland. We found qualitatively similar responses of host taxa (Figure 6) to the same type of climatic variables, which were found to affect host use (Figure 3) in the large, pan‐arctic data: warming summer time temperatures translate into increases in lepidopteran herbivory and decreases in the abundance of dipteran pollinators (Figure 1). Warmer autumn soil has the opposite effect of decreasing lepidopteran herbivory and increasing muscid fly abundance. Together these factors result in opposing temporal population patterns in herbivore and pollinator populations (Figure 6a). For dipteran pollinators, and specifically for the muscid flies, part of the effects here attributed to soil temperatures is perhaps more likely reflective of other factors covarying with temperature, such as the amount of snow and more specifically soil water content during summer. Hence, soil temperature serves as a proxy of other corollaries of climate change. A more detailed analysis of the muscid fly communities at Zackenberg was provided by Loboda et al. (2018) and reviewed in Gillespie et al. (2019). Their species‐specific assessment showed that some species fail to recover after bad years, thus altering the community composition. Such dynamics could explain why our family‐level analysis underestimates fly abundances in the early years and overestimates them for the most recent years, and highlights the need for taxonomic resolution in examining community responses to environmental change. Reponses to soil water may also account for an added and somewhat counterintuitive pattern detected: that a higher number of trapping days were associated with fewer fly individuals caught per season. Here, the underlying reason is simple: long trapping seasons are generated by early snow melt. Thus, the relationship between time period and fly counts is not a causal one, but generated by an indirect association between the length of the local trapping season, the habitat, and the fly abundances caught.

Climate change effects on large spatial scales

The key insights delivered here build on a new type of space‐for‐time approach, where regions are characterized not only by different climates but also by the rate and mode of recent climate change, that is, a space‐for‐change approach. We further support these inferences by comparison of our spatial data with local time series (Figure 2). In support of our overall interpretation that recent climate change has affected the community composition of parasitoids, we stress two considerations emanating from this unique combination of data sources (Figure 2). First, the rate of temperature change during summer and winter explains more of the observed variation than does any metric of mean climatic conditions examined. Second, the time series data (Figure 2) corroborate the trends observed for two different host arthropod guilds (herbivores and pollinators) and produce a strikingly similar pattern of correlations with environmental variation.

Overall, the possibility of measuring climate change consequences by comparing regions of contrasting recent climate changes histories offers hope for filling crucial knowledge gaps. Given the rapid changes in arctic climates, there is an urgent need for functionally meaningful descriptions of parasitoid communities in the Arctic, and of arctic arthropod communities in general (Gillespie et al., 2020). Systematically collected long‐term datasets of insect abundance from the Arctic are rare, limiting current insights into how climate change affects community structure, dynamics, and functioning. To our knowledge, temporal change in arctic parasitoid community composition over time has been assessed in only two instances: by Fernandez‐Triana et al. (2011) at Churchill, Canada and by Timms et al. (2013) on Ellesmere Island, Canada. The two studies found contrasting results, with the study by Fernandez‐Triana et al. (2011) reporting high species turnover at a subarctic site, whereas Timms et al. (2013) observed practically no change in species composition on an isolated island in the High Arctic. This pattern is in part expected also based on the sites' climate change histories, which Churchill having experienced faster change. However, these studies both made comparisons across datasets resulting from different methodologies, limiting their analysis to a presence/absence approach. The functional consequences of these two different community level responses remain unknown, but we believe that future studies utilizing the space‐for‐time approach presented here may help resolve them.

Both mean climate and recent change shape functional community composition

In our exploration of parasitoid community structure, we targeted two different aspects of community composition among the parasitoids of Lepidoptera: the proportion of idiobiont species, which is likely to change more slowly over time, and the proportion of idiobiont individuals, which can change at a timescale of only some year(s) (Figure 1). We show that against the backdrop set by regional climate, the idiobiont strategy is more common than expected in rapidly warming areas. We interpret this as a sign of recent community change (question 3). While ideally we would have validated this finding against changes over time in the Zackenberg time series, this remains out of reach—as the Zackenberg samples remain to be identified to a taxonomic level allowing assignment of taxa to idiobionts versus koinobionts (but see Ji et al., 2019).

In terms of responses at species level, differential responses among idiobionts and koinobionts are consistent with what we know about their biology. As such, the number of idiobiont species drops more sharply with winter temperature than does the number of koinobiont species (Quicke, 2012; Timms et al., 2016). This pattern fits well with the trade‐offs association with parasitoid life‐history strategies, where some koinobionts benefit from the sheltered overwintering burrows of their hosts. This effect was found in taxa parasitizing predominantly Lepidoptera, whereas such trends were not detected in the dominance of idiobiont species within the species pool of Diptera‐using parasitoids. Such a contrast makes sense, as dipteran larvae tend to relatively less mobile, and will hardly spend much time seeking for suitable sites for pupation. Thus, koinobiont and idiobiont parasitoids of Diptera may be faced with approximately the same overwintering conditions. The strong negative response to winter warming observed in the proportion of idiobiont individuals among parasitoids of Diptera is likely explainable by stronger associations between parasitism strategy and different host lifestyles. As an example, the most abundant parasitoids of fungus gnats (Diptera: Mycetophilidae) are koinobiont. Thus, among parasitoids of Diptera, host abundance changes are liable to override the independent effect of parasitism strategy, unlike in Lepidoptera‐using parasitoids, in which koino‐ and idiobiont species commonly share host species (in our dataset).

In terms of responses at the population level among the parasitoids of Lepidoptera, winter temperature was a key driver restricting the relative abundance of idiobionts. However, recent changes in climatic conditions emerged as a potent force overriding mean conditions. The more a location has warmed during the past 18 years, the higher the fraction of idiobiont lepidopteran parasitoids—relative to what would be expected based on the regional mean climate. The trend is likely explained by species traits associated with koino‐ or idiobiontism. Not only will warming winters favor idiobiont survival (Hance et al., 2007) but also an increase in the summer temperatures may improve their host searching (Gu & Dorn, 2001; Menon, Flinn, & Dover, 2002). In contrast, koinobionts are more likely to suffer from phenological mismatches with their hosts as the time windows conducive to development get shorter (Duan, Jennings, Williams, & Larson, 2014; Van Nouhuys & Lei, 2004). Warming may also disrupt host immunosuppression, on which species with a koinobiont lifestyle rely (Seehausen et al., 2016; Seehausen, Régnière, Martel, & Smith, 2017). This would lower the overall mortality to parasitism, but also disproportionately affect koinobionts. All in all, we consider the idiobiont‐to‐koinobiont ratio to be an accurate, temperature‐sensitive metric for tracking functional community change among parasitoid communities.

Determinants of herbivory in the Arctic

Judging from both the large‐scale data and the especially the local Zackenberg time series, our results point to increasing arctic herbivory with warming summers (question 4) and the opposite with warming winters. These findings are consistent with those from a pan‐arctic survey of background leaf herbivory, which revealed summer temperature as a key driver of insect herbivory levels across the whole arctic plant community (Rheubottom et al., 2019) although with large variability. This pattern was previously observed in more detailed study on dwarf birches (Betula glandulosa‐nana complex; Barrio et al., 2017). Contrasting with our results, the largest effects in that study were found at locations with warmest mean temperatures. This difference may perhaps be due to contrasting features in the plant–herbivore systems examined. For the dwarf birch system, a positive correlation between snow cover and levels of herbivory has been suggested (Torp, Olofsson, Witzell, & Baxter, 2010). In contrast, for Dryas, herbivores seem to benefit from snow‐free conditions, either directly due to lower energy consumption during overwintering (Bowden et al., 2015) or indirectly due to a prolonged growing season. Furthermore, the latitudinal abundance patterns of Dryas and Betula are contrasting: as, with Dryas increasing northwards, and Betula declining. Finally, the primary herbivores of Dryas are Lepidoptera (which presumably decline in number northward), whereas dwarf birches host various herbivore taxa, some of which remain relatively diverse and abundant at higher latitudes (i.e., sawflies; Kouki et al., 1994). In conclusion, climate change may have different consequences for different plant–herbivore systems, depending on the relative effects of temperature, snow, and precipitation. This complexity may explain why effects of summer warming were not as clear across the Arctic as they were in a single locality experiencing particular type of change.

We reveal a joint impact of mean climate and recent climate change, with important repercussions for arctic ecosystem functioning. Consistent with the general pattern of a community shift from dipteran to lepidopteran host use with a change in climate, we show that the same climatic factors that benefit lepidopteran herbivores are detrimental to dipteran pollinators. As our results are mostly derived from relatively undisturbed arctic sites, they offer a rare view of climate‐driven processes largely untainted by other (anthropogenic) impacts. As such, they provide insights into climate‐driven effects on biotic interactions, and point to a general instability of community structure and to increasing trends in herbivore damage. We thus call for further arthropod community studies in the Arctic, to better understand community‐level impacts at the leading edge of climate change. In addition, we highlight the benefit of using distributed study designs drawing on a network of collaborators, thereby allowing standardized sampling across diverse field sites.

Supplementary Material

Click here for additional data file.

Supplementary Material

Click here for additional data file.