Accounting for species interactions is necessary for predicting how arctic arthropod communities respond to climate change

1 –––––––––––––––––––––––––––––––––––––––– © 2021 The Authors. Ecography published by John Wiley & Sons Ltd on behalf of Nordic Society Oikos This is an open access article under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited. Subject Editor: Eric Post Editor-in-Chief: Miguel Araújo Accepted 19 February 2021 44: 1–12, 2021 doi: 10.1111/ecog.05547 44 1–12


Introduction
A central goal of current biodiversity research is to better understand and predict biodiversity response to climate change. Most research on how climate change will affect future ecological communities has focused on the link between the abiotic changes caused by climate warming and subsequent changes in species distributions and abundances (Pereira et al. 2010, Pacifici et al. 2015. Species extinctions and changes in species' distribution and abundance are currently modifying the composition and functioning of ecological communities (Parmesan 2006, Chen et al. 2011, Blowes et al. 2019. However, the ongoing changes are modulated not only by the direct impacts of changing climatic conditions on individual species but also indirectly, through Accounting for species interactions is necessary for predicting how arctic arthropod communities respond to climate change cascading effects of species interactions (Petchey et al. 1999, Memmott et al. 2007, Tylianakis et al. 2008, Harley 2011. Cascade effects follow when the abundance of a species in a local community changes in response to e.g. climate warming, influencing the occurrences or abundances of other species it interacts with. Therefore, the joint assessment of climate and interspecific interactions is of fundamental importance for improving our current understanding and our predictions of how climate change affects ecological communities (Gilman et al. 2010, Van der Putten et al. 2010, Blois et al. 2013, Scheffers et al. 2016.
In the Arctic, the rate of temperature increase is nearly double the global average (Kattsov et al. 2005, IPCC 2014, Post et al. 2019. This is dramatically changing abiotic conditions, with important impacts on phenology, especially at the onset of the growing season (Høye et al. 2007, Wheeler et al. 2015, Kankaanpää et al. 2018. As a consequence, plants and animals adapted to longer growing seasons as well as species originating from lower latitudes, or from warmer locations nearby, are becoming more abundant (Sturm et al. 2001, Post et al. 2009, Elmendorf et al. 2012, Elmhagen et al. 2015. Conversely, two types of species are predicted to decline in particular: species adapted to the short but intense growing season of the Arctic, and species sliding out of synchrony with their key resources Forchhammer 2008, Hegland et al. 2009). Related to the latter, the effects of climate warming on arctic communities are intensified by their impacts on species interactions. Changes in phenology can decouple interacting species, causing 'phenological mismatch' (Renner and Zohner 2018), and such mismatches can have serious consequences, including population declines of herbivores Forchhammer 2008, Ross et al. 2017), specialist predators (Gilg et al. 2009, Schmidt et al. 2012a) and pollinators (Høye et al. 2013, Schmidt et al. 2017, or -mutatis mutandis -the escape of herbivores from enemies and increased plant damage (Kankaanpää et al. 2020(Kankaanpää et al. , 2021. Arctic terrestrial interaction webs are dominated by arthropods (Wirta et al. 2015(Wirta et al. , 2016. As recently uncovered by molecular tools, the diversity of arctic arthropods exceeds what was previously assumed (Wirta et al. 2015(Wirta et al. , 2016. The short lifespans, strong functional responses to temperature and narrow phenological niches of arthropods make them particularly responsive to warming (Deutsch et al. 2008). Indeed, the few studies of temporal trends in population abundance already suggest that arctic arthropod communities are experiencing severe changes (Gillespie et al. 2019, Kankaanpää et al. 2020. Among soil and foliage-dwelling arctic arthropods, herbivores (Hemiptera and Lepidoptera) and parasitoids (Hymenoptera) have increased locally in abundance, whereas detritivores (Collembola and Acari) have decreased (Koltz et al. 2018a). While spider communities (Araneae) seem not to have experienced drastic changes overall, individual species have declined in abundance (Bowden et al. 2018). Among flower-visiting arctic arthropods, muscid flies (Diptera) have drastically decreased in abundance (Loboda et al. 2018). Also the abundances of aquatic arthropods such as chironomids (Diptera) are changing, with abundances predicted to increase with rising temperatures (Engels et al. 2020, but see Høye et al. 2013).
Some of the most frequent trophic interactions among arctic arthropods involve spiders feeding on herbivorous arthropods (Wirta et al. 2015, Eitzinger et al. 2019, and parasitoid wasps and flies feeding on lepidopteran and dipteran larvae (Várkonyi andRoslin 2013, Kankaanpää et al. 2020). Given these strong trophic interactions, climateinduced change in arthropod diversity and abundance can be expected to cascade through these high-latitude communities (Koltz et al. 2018b, but see Visakorpi et al. 2015). Furthermore, changes in arthropod communities can have drastic consequences on ecosystem functioning (Schmidt et al. 2017). For instance, a large proportion of tundra plants is insect-pollinated (Kevan 1972), with muscid flies identified as the main pollinators of some dominant plant species (Gillespie et al. 2016, Tiusanen et al. 2016. Also, while plant consumption by arctic arthropods is rather low , it is still of the same magnitude as herbivory by large mammals (Mosbacher et al. 2016). Thus, predicted outbreaks of herbivore arthropods in the Arctic can have important consequences for plant consumption (Jepsen et al. 2008) and ecosystem primary productivity (Lund et al. 2017). Finally, arthropods are the main food resource for many migratory birds during their breeding season, and changes in arthropod availability may thus affect arctic bird reproduction and growth (Meltofte et al. 2007, van Gils et al. 2016. Despite this evidence for pervasive trophic interactions among arctic arthropods, coupled with the high sensitivity of arthropod population dynamics to climate change, the role of trophic interactions to their responses to global climate change is still largely unknown. Here we aim to fill this knowledge gap by asking whether we can use information on species interactions to improve our predictions of arctic-arthropod population dynamics. For this purpose, we use a 14-year-long time series of arthropod communities from Zackenberg, Northeast Greenland collected by BioBasis (Schmidt et al. 2012b) as part of the Greenland Ecosystem Monitoring (GEM) program (<https://data.g-em.dk/>, accessed 9 Aug 2020). The time series consists of physical samples of all trophic levels within the local arthropod community. While most specimens have previously been classified to higher taxonomic ranks only (typically families, with species-level identifications available for only a minority of the specimens; Schmidt et al. 2016 for exact categories), in this paper we achieve species-level resolution for all samples by applying a mitogenomic-mapping pipeline that we recently developed (Ji et al. 2020). To elucidate the impact of species interactions on the population dynamics of arthropods, we compared the predictive performances of four alternative joint species distribution models, which differed in whether and how interactions were taken into account (Fig. 1).

Study system and sampling methods
The data were collected as part of the GEM program in Zackenberg, a High Arctic site located in Northeast Greenland (74°28′ N, 20°34′ W). Zackenberg is characterized by a continental climate with mean monthly temperatures ranging from −20 to +7°C and an annual precipitation of 260 mm. The vegetation consists of tundra species, of which the arctic willow Salix arctica, the arctic bell-heather Cassiope tetragona and the mountain avens Dryas spp. are some of the most abundant. The growing season has lengthened by ca 50% during the study period, from a snow-free period lasting from end-June to end-August in the late 1990s to a period lasting from mid-June to the beginning of September in the 2010s (Kankaanpää et al. 2018).
Arthropods were sampled weekly during the growing season from 1997 to 2013, using yellow-pitfall traps located in a mesic heath habitat. While the original number of traps monitored by the BioBasis program is larger, we focused on traps from one of the six still-operational arthropod monitoring sites within the larger landscape (Ji et al. 2020). Specifically, the material derives from three of four traps (A-C) collected at trapping station ART3 (for maps and details, Schmidt et al. 2012b), the traps being located 5 m from each other. Unfortunately, the 2010 samples were lost in transit from Greenland, so data from this year are missing.

Sample processing, sequencing and mitogenomic mapping
In a previous methods publication (Ji et al. 2020), we described the SPIKEPIPE mitogenome-mapping pipeline, which converts arthropod bulk samples into data on species presence-absence and DNA-abundance. Below, we briefly summarize the steps performed by Ji et al. (2020), followed by a description of how we supplemented the data matrices generated by Ji et al. (2020) with information from other sources. More detailed descriptions of sample processing, sequencing and bioinformatics are in the Supporting information and Ji et al. (2020).
The arthropod samples were preserved in ethanol at room temperature in the Museum of Natural History, Aarhus, Denmark. Over the years, the samples have been sorted to subsamples of higher taxonomic levels (typically families). In 2016 and 2017, we non-destructively extracted genomic DNA from these subsamples following a modified procedure from Gilbert et al. (2007), pooled them into the original trap-week samples, added a fixed-aliquot DNA-spike as a standard, shotgun-sequenced the DNA from each sample, and mapped the output reads to a mitogenome reference database. Each trap-week sample was individually libraryprepped and sequenced at the Earlham Institute in Norwich, UK. The mitogenome reference library was constructed from the voucher collection established by Wirta et al. (2016). We used the DNA-spike in each sample, technical replicates across sequencing runs, and a mock-sample-estimated Figure 1. Conceptual illustration of the four time-series models compared in this study. Left. Species data are included in the statistical models as presence-absence of species in a particular trap-week (w) of a particular year (Y) (when used as a response variable, top), and as abundance measured by the proportion of trap-weeks during which the species was recorded in a previous year or absolute abundance measured by DNA (when used as a predictor, bottom). Data on climatic variables (mean air temperature and snow depth) are included as predictors (top). Right. The models of increasing complexity assume that the temporal dynamics of species communities are determined by an increasing number of factors: 1) abiotic environmental variables (sampling week temperature and previous summer temperature or previous winter snow depth) only (abiotic model), 2) abiotic environmental variables + intra-specific interactions (density dependent model), and abiotic environmental variables + intra-specific interactions + inter-specific interactions, the latter structured either 3) at the species level (species-interactions model) or 4) at the trophic-group level (trophic-interactions model). For simplicity, in this figure the models are illustrated from the point of view of predicting the dynamics of one focal species, whereas in the joint species distribution models, all species in the community data are simultaneously focal species. threshold for minimum-acceptable mapping-coverage to achieve accurate data on species presence-absence as well as DNA-abundance.

Community data
The species data matrix derived from mitogenome mapping was complemented with data from two other sources. First, samples from the earliest and latest parts of the summer often comprised only a few individuals. For these, we DNAbarcoded each individual with PCR primers LCO1490 and HCO2198 (Folmer et al. 1994), identified them to species using the BOLD database (Ratnasingham and Hebert 2007), and added their occurrences to our data matrix. These DNAbarcode data include only the species also detected by mitogenome mapping. Second, we added data on Diptera species from subsamples that had been identified in an independent morphological study (Loboda et al. 2018) but had not yet been returned to the collection, and which thus could not be sequenced. This added three new Diptera species to the data, as well as new occurrences of species already detected by the mitogenomic sequencing. Altogether, the community data compiled consisted of 81 arthropod species identified in 542 trap-week samples. Given that this dataset represents arthropods co-occurring in space (trap) and time (week), we consider the species found in the trap-week samples to belong to the same community (Fauth et al. 1996).
To include functional information in the form of the trophic levels represented in the community data, each species was classified by its non-adult feeding type: herbivore (9 species), omnivore (1 species), parasitoid (18 species), predator (25 species) or species that were either saprophages or detritivores (28 species). These classifications are provided in the Supporting information.

Climate data
From GEM's ClimateBasis program, we downloaded average hourly soil and air temperatures (the latter measured 2 m aboveground), as well as three-hourly data on snow depth, spanning the study period. As the soil and air temperatures were highly correlated (Pearson correlation r = 0.90), only airtemperature data were used. We calculated from these data three climatic predictors: 1) the mean temperature during the week at which the arthropod sampling event was conducted (henceforth called sampling week temperature); 2) the mean temperature during the previous summer (averaged over weeks 23-32; henceforth called previous summer temperature); and 3) the mean snow depth during the previous winter (averaged over weeks 33-52 and 1-22; henceforth called previous winter snow depth). The sampling week temperature predictor was aimed to capture insect activity and thus detectability (Høye and Forchhammer 2008), whereas the previous-year summer temperature and previous-year winter snow depth predictors were aimed to capture the influence of climatic factors on population growth rate. Previous summer temperature has been found to affect arthropod population growth rates (Koltz et al. 2018a), while snow depth is associated with the timing of snow melt, which affects plant phenology and consequently arthropod population dynamics (Høye et al. 2013).

Analyses of species richness
As a measure of species richness, we used the yearly averages of the number of species identified per trap-week sample. We used linear regression to quantify trend in species richness a function of year. To examine whether the trend differed among trophic groups, we repeated the analysis separately for each trophic group. We further examined if there has been a systematic trend in temperature by regressing temperature against the linear effect of year. To account for possible temporal autocorrelation in the model residuals, we supplemented the baseline linear model analyses with an ARIMA (autoregressive integrated moving average) model. We implemented the linear model with the 'lm' function in R (<www.r-project. org>) and assessed the significance of the linear trend by the p-value. We implemented the ARIMA model with the 'arima' function in R and assessed the significance of the linear trend by the comparison of AIC values between model variants that included and excluded the linear trend of the year.

Analyses of population dynamics
To identify the drivers of population dynamics, we applied a joint species distribution model (Warton et al. 2015) via the hierarchical modelling of species communities (HMSC) (Ovaskainen et al. 2017b, Ovaskainen and Abrego 2020), using a time-series approach similar to that of Ovaskainen et al. (2017a). While the data were generated at a weekly resolution, our main interest was in the yearly transitions, as most arctic arthropods have an annual or multiannual life cycle, and thus the population dynamics of most of the arthropod species in the dataset take place on an annual basis. Thus, the predictor variables of previous-year species abundances and previous year climatic conditions were defined as yearspecific averages (Fig. 1). The response variable was however measured at the resolution of species-trap-week, as this allowed us to account for both phenological variation and the influence of weather conditions on insect activity and thus detectability (Høye and Forchhammer 2008).
Concerning species abundances in the previous year, we considered two different predictors: 'prevalence-based abundance' and 'DNA-based abundance' (Fig. 1). Prevalencebased abundance was defined as the fraction of trap-week sampling units occupied by a given species over the previous year. We considered this measure to be an ecologically relevant proxy of species abundance, because if a species is consistently present in many samples from a given year, then this species is likely to be more available for species interactions. While the DNA abundance is in theory more informative than merely prevalence, it includes species-specific biases (e.g. mitochondrial copy number, individual biomass), for which reason we scaled the estimates of Ji et al. (2020) to zero mean and unit variance within each species. The DNA abundance can also be expected to be more prone to random variation due to e.g. outliers generated by an exceptionally large DNA yield. Further, for part of the data, DNA abundance could not be recovered (i.e. for Sanger sequenced and non-sequenced individuals), so these samples represent missing data for DNA abundance.
We fitted a probit-model where the response variable was the vector of presences and absences of all species in a particular trap in a particular week in a particular year. In all models, we included as predictors the linear and second order effects of the week (to model species phenology), and the sampling week average temperature to capture variation in thermal conditions that influence insect activity and thus detectability (Høye and Forchhammer 2008). Furthermore, we included either previous-year summer temperature or previous-year winter snow depth as the climatic predictor influencing insect growth rate, either by directly affecting insect thermal tolerances and/or vegetation phenology (Høye et al. 2013, Nabe-Nielsen et al. 2017, Kankaanpää et al. 2018, Koltz et al. 2018a). The reason for fitting alternative models rather than including both predictors into the same model was that the time series contained in total 14 yearly transitions only, and thus it was not possible to include a large-number of year-specific predictors. The abiotic model included no additional predictors. In the density-dependent model, we added as a predictor the abundance of only the focal species in the previous year (Fig. 1). In the trophic-interactions model, we likewise included the abundance of the focal species in the previous year, as well as four additional predictors measuring the abundances of herbivores, parasitoids, predators and saprophages and detritivores in the previous year, each averaged over the species within each trophic group. In the species-interactions model, we included as predictors the species-specific abundances of occupied trap-weeks in the previous year for all species, thus modelling both within-and among-species density dependence at the species level.
To account for the spatial and temporal non-independence of the sampling units, we included trap (three levels of trap A, B and C) and year (14 levels) as random effects in all models. To yield community-level inference, we further added a hierarchical model structure indicating the trophic group of the species. In the HMSC framework, traits are used as predictors for how the species respond to fixed effects (Ovaskainen et al. 2017b).
In the HMSC analyses, we included only the 52 species that were present in at least 10 trap-week samples (Supporting information). Furthermore, while the entire data consisted of 542 trap-week samples, in the HMSC analyses we included only 467 samples, as we needed to exclude data from the year 2011, for which we could not compute the biotic predictors since the data from 2010 were missing. We fitted the models using the R-package Hmsc 3.0 (Tikhonov et al. 2020).
In the species interaction model we specifically avoided over-parameterization by applying the so called 'sparse interaction model' of Ovaskainen et al. (2017a). In this approach, the level of variable selection is tuned by choosing the prior probability by which each variable (in this case the abundance of a specific species in the previous year) is included in the model, i.e. by which the regression coefficient is different from zero. To test the sensitivity of the results with respect to this choice, we varied the prior probability from 0.01 to 0.1, 0.5 and 1.0. We then chose the model with the highest predictive performance, i.e. the one that was best able to capture the signal while avoiding overfitting to the noise. We evaluated the level of over-parameterization by comparing the explanatory power (predictions based on models fitted to all data) with the predictive power (predictions based on a cross-validation approach. For technical details about model fitting and model comparison, see Supporting information.

Results
Over the study period, average summer temperature at Zackenberg increased by 2.0°C (from 1.28°C to 3.28°C, p = 0.012 in the linear model and deltaAIC = 5.4 in the ARIMA model, Fig. 2a), and arthropod species richness drastically decreased by 50% (from a yearly mean of 8.9 species to 4.4 species per trap-week, p = 0.003 for the linear trend, Fig. 2b). The overall species richness increased to unusually high level in 1999, then crashed during 2000-2001, after which it first partially recovered and then continued decreasing until the end of the study period (Fig. 2b). Changes in species richness were group-specific, indicating that the food-web structure changed over time (Fig. 2c-f ). Predators and saprophages/detritivores decreased the most steeply over the study period, whereas for herbivores and parasitoids, the dynamics show a fluctuating pattern, with no overall trend. Since most parasitoids are in the family Ichneumonidae, most predators in Muscidae and most detritivores in Chironomidae (Supporting information), trends in these families match trends at the level of their corresponding trophic-group (Supporting information).
Weekly variation in presence-absence was satisfactorily predicted by all models, as shown by the mean AUC value over the species ranging from 0.82 to 0.88 for all model variants. However, our focus was at a yearly resolution, i.e. at assessing how well the models were able to predict in which years species abundances were high or low. At this level, the explanatory powers of the models, measured as the correlation between predicted and observed yearly abundances, ranged from 0.61 to 0.99 (Table 1). As may be expected, the explanatory powers generally increased with increasing model complexity, since more complex models have more parameters that can be estimated. To account for this potential effect of overfitting, we evaluated the predictive powers of the models through a cross-validation procedure where we masked the data from each year for which a prediction was made. Compared to the explanatory powers, the predictive powers were much lower (Table 1), highlighting the difficulty of ecological prediction. The predictive powers of the abiotic models were even negative, reflecting the fact that when excluding a focal year from the model fitting, the predictions on populations abundances may be opposite to reality (i.e. predicted to be high in a year that were low in reality, and vice versa).
When we added interspecific interactions, either at the level of pairs of trophic groups or pairs of species, predictive power turned from negative to positive (Table 1). The highest predictive powers were achieved by the pairwise species-interaction model, which estimates the full pairwise species-to-species interaction matrix of the yearly transitions. In terms of how species abundance in the previous year was measured, species prevalence turned out to be a better predictor for cross-validation than DNA abundance for all model variants, for which reason we present below results only for the model using prevalence-based abundance.  Table 1. Explanatory and predictive powers of the four joint species distribution models varying in their abiotic predictors (previous summer temperature and previous winter snow depth). We measured the yearly explanatory power (r) by calculating the Spearman rank correlation between the predicted and observed arthropod occurrences, both averaged into yearly prevalences. The yearly predictive power (CV-r) was measured by year-based cross-validation, i.e. by assessing how well the models fitted to data excluding each focal year were able to predict the arthropod communities. Note that all the models include also the explanatory variables of sampling week temperature and effect of the week. We measured the weekly explanatory power (AUC) by calculating the area under the curve between the predicted and observed arthropod occurrences. The weekly predictive power (CV-AUC) was measured by 10-fold cross-validation. The results shown here are based on prevalence-based abundance as the predictor, and the results for the models using DNA-based abundance are given in the Supporting information.

Model
Previous

891
The results for DNA-abundance are shown in the Supporting information. For the abiotic model, the model variant including previous summer temperature as the climatic predictor performed better (less negative predictive power) than the variant including previous winter snow depth (Table 1). Interestingly, the model which achieved highest predictive power was the species interaction model which included the previous winter snow depth as climatic predictor. For the model variant including previous summer temperature as the predictor, accounting for species interactions increased the predictive power by 0.18 (from −0.08 to 0.1). For the model variant including previous winter snow depth, accounting for species interactions increased the predictive power much more, by 0.36 (from −0.20 to 0.16). This result suggests that species abundances in the previous years and snow depth in the previous winter bring complementary information.
As expected, the difference between predictive and explanatory power was greatest for the most complex models, as these have more parameters and thus a higher risk for overfitting. However, in the species interaction models the predictive power did not improve when applying the so called sparse interactions model (Supporting information). Hence, the results shown here correspond to the model without variable selection.
When assessing how the predictive powers varied among species, we observed that in the model variant including previous summer temperature, the predators achieved generally the highest and the herbivores the lowest predictive power (Supporting information), whereas no systematic differences were observed among the trophic groups in the model variant including previous winter snow depth as the climatic predictor (Supporting information). In both model variants, the predictive power increased with species prevalence in the abiotic model but not in the species interactions model (Supporting information), suggesting that accounting for species interactions increased the predictive power especially for the rare species.
Variance partitioning among the explanatory variables showed that arctic arthropod communities are characterized by marked seasonal variation (week and its square explained more than half of all explained variation in all but species interaction models, Fig. 3). Climatic conditions also affected arthropod community dynamics, but only to a lesser extent (sampling week temperature and previous summer temperature/previous winter snow depth explained altogether 12-17% of the variation in all but species interaction models). Concerning the biotic predictors, density-dependence within species explained only a minor part (3%), trophic interactions a substantial part (20-22%) and species interactions a major part (51%) of the variation.
For most species, the responses of arthropod species to the previous week's temperature were positive with at least 90% posterior probability, reflecting the positive effect of withinseason temperature on arthropod activity and thus detectability (Fig. 3). In all four models, the response of arthropod species was negative to the second-order effect of week, reflecting the peak in arthropod's occurrence at intermediate weeks. Within-species density dependence was predominantly positive in the few cases that showed an effect (Fig. 3). In the trophic-interactions model, parasitoid and predator abundances positively influenced many species, especially saprophages and detritivores (Fig. 3). In contrast, saprophage and detritivore abundance influenced many species negatively, especially herbivores (Fig. 3). Herbivore abundance influenced parasitoids positively and saprophages and detritivores negatively (Fig. 3). In the species-interactions model, no interactions were supported with over 90% posterior probability, demonstrating the difficulty of inferring exactly which species pairs are driving the community-level trends.

Discussion
Developing accurate predictive models of biodiversity change is a priority for global-change science. A topical issue in formulating such models is whether and how species interactions should be accounted for, especially when modelling species-rich communities (Gilman et al. 2010, Blois et al. 2013). Our study of arctic arthropods demonstrates that accounting for interspecific interactions either at the species or trophic group level indeed improves the predictive performance of time-series models.
In demonstrating the importance of biotic interactions for understanding community dynamics, we derive support from and build on several studies that have improved predictive performance by including the abundances of one or a few interacting species as predictors in single-species distribution models (Araújo and Luoto 2007, Heikkinen et al. 2007, Wisz et al. 2013, le Roux et al. 2014, Mod et al. 2015. In our study, we have extended this approach to a communitywide perspective, showing how to account for a complex interactive network via joint species distribution models that describe community-level dynamics. Interestingly, accounting for within-species density dependence did not improve predictive power, whereas accounting for inter-specific interactions improved the predictive powers of the models. Given the imperative need to scale up from species-level to community-level projections of biodiversity change, these results represent a major step forward for predictive ecology. Arthropods are fundamentally important ecosystem components, both in species numbers and in contributions to ecosystem functioning (Schmidt et al. 2017). We report an alarming rate of arthropod decline in a High Arctic site, particularly in predators and saprophages/detritivores. Our study site is essentially untouched by human activity, except for climate warming, which is progressing twice as quickly in the High Arctic as the global average (Kattsov et al. 2005, IPCC 2014, Post et al. 2019. Declines in arctic arthropod species richness and abundances of several taxonomical and functional groups have previously been reported in relation to climate warming (Høye et al. 2013, Koltz et al. 2018a, Loboda et al. 2018, Gillespie et al. 2019. This is the first study presenting comprehensive community-level results with species-resolution data, allowing us to compare the rates of decline among groups. Our observation of declining trends among arthropods echoes those of Koltz et al. (2018a), who partially used samples from the same study area. A key difference is in the information inferred from species-versus group-level data. Where Koltz et al. (2018a) used aggregate counts of higher taxonomic ranks, we used data on each species in the community to build support for emergent change.  That is, rather than assuming that a group-level response reflect consistent change among the species in the group, we directly evaluated this assumption by observing species-level changes within and among groups. Taken together, the results support the interpretation that the food-web structure of arctic arthropod communities has changed markedly over the last decades (Kankaanpää et al. 2020). From our results, it is clear that the trophic-web structure is radically changing in the High Arctic, and that this affects arthropod community composition much more profoundly than if impacts derived solely from the direct effects of climatic changes. Especially predators and saprophage/ detritivores have declined disproportionally compared to the other groups (this study, Koltz et al. 2018a). Our result show that the trajectory of community change will be driven not by each species responding individually to the changing environment, but by complexed lagged effects between trophic groups. As such, our findings add evidence to reports on insect decline (Cardoso and Leather 2019, Sánchez-Bayo and Wyckhuys 2019) and add to concerns regarding subsequent changes in ecosystem functioning. Working in Alaska, Koltz et al. (2018b) found changes in climate to reverse top-down effects of predators on belowground ecosystem function, mediated by changes in mesofauna (springtails [Collembola] and mites [Acari]). Our study shows disproportionate changes among macrofauna with a detritivorous larval stage (Diptera, Insecta), suggesting an Arctic-wide pervasive impact of climate change on detritivores.
The result that different arthropod trophic levels are differently affected by climate warming in the Arctic is consistent with previous reports of shifts between specific trophic levels. With respect to plant-herbivore interactions, experimental (Roy et al. 2004, Liu et al. 2011, de Sassi and Tylianakis 2012, Birkemoe et al. 2016) and observational studies (Barrio et al. 2017, Rheubottom et al. 2019, Kankaanpää et al. 2020 report changes in arthropod-driven herbivory. With respect to predator-detritivore interactions, experimental studies suggest that effects of climate change can cascade through other trophic levels, thereby altering critical ecosystem functions (Koltz et al. 2018c). In our study system, we note that a previous attempt to detect trophic cascades from data on densities of a single predator spider species Pardosa glacialis was unsuccessful, a finding attributed to the complexity of the local food web and the many outflows from an increased predation pressure (Visakorpi et al. 2015). Likewise, the current result that accounting for density dependence did not improve model predictions may reflect the fact that capturing density dependence in arthropods is challenging, due to highly variable population dynamics (Hanski et al. 1990) and/or to the fact that trapping success is influenced both by abundance and activity. In this study, we included a higher level of complexity in the analyses by jointly analysing 52 arthropod species, classified by trophic group, enabling us to detect trophic-cascade effects across the arthropod community.
Some of the trophic-group effects are easy to interpret, such as the positive influence of the previous year's herbivore abundance on parasitoids. Other results are more difficult to interpret, such as the positive influence of parasitoids and predators on saprophages/detritivores, and the negative influence of saprophages/detritivores on herbivores. Whether these effects are due to direct or unknown species interactions or some (delayed) responses of the species to environmental conditions not included in our modelling is an important question for future studies. Important but unaccounted abiotic predictors may for example include snow and microclimatic conditions (Kankaanpää et al. 2018(Kankaanpää et al. , 2020a. Thus, we note that variation assigned to trophic groups may in some cases represent covariation with unmeasured environmental covariates rather than the direct effects of biotic interactions.
However, we stress that identifying the exact causality is not necessary for our core inference: that accounting for interspecific associations is necessary for successfully predicting how arthropods respond to climate change. Regardless of whether the parameters of the species and trophic interaction models relate causally to biotic interactions or to synchronous responses to unmeasured environmental conditions, they were successful in improving predictive power. However, this does not imply that identifying the causal components is not important. In fact, a model that makes better predictions 'for the wrong reasons' will be able to successfully predict future patterns only for as long as the relationship between the causal factors and the apparent patterns of species interactions also remain constant into the future.
We found that models using only climatic predictors generated predictions that were opposite to observations. This adds an important caveat to the many studies using correlative, climate-envelope approaches to model the response of biodiversity to climate change (Pacifici et al. 2015, Warren et al. 2018, Trisos et al. 2020. Our results also illustrate the absolute difficulty of predicting species communities, given poorer predictive performance versus explanatory performance, and attest to the increased risk of overfitting when model complexity increases with the inclusion of more species and environmental predictors. While our conclusions are based on arctic arthropods, since climate warming and species interactions are ubiquitous, we suggest that our approach serves as a conceptual and methodological template for deriving predictive models for other organism groups and ecosystems.

Data accessibility statement
The data and code for replicating the results are deposited in Dryad Digital Repository <https://doi.org/10.5061/dryad. cc2fqz65p>.