Temperature effects on the temporal dynamics of a subarctic invertebrate community

1. Climate warming is predicted to have major impacts on the structure of terrestrial communities, particularly in high latitude ecosystems where growing seasons are short. Higher temperatures may dampen seasonal dynamics in community composition as a consequence of earlier snowmelt, with potentially cascading effects across all levels of biological organisation. 2. Here, we examined changes in community assembly and structure along a natural soil temperature gradient in the Hengill geothermal valley, Iceland, during the summer of 2015. Sample collection over several time points within a season allowed us to assess whether temperature alters temporal variance in terrestrial communities and compositional turnover. 3. We found that seasonal fluctuations in species richness, diversity, and evenness were dampened as soil temperature increased, whereas invertebrate biomass varied more. Body mass was found to be a good predictor of species occurrence, with smaller species found at higher soil temperatures and emerging earlier in the season. 4. Our results provide more in-depth understanding of the temporal nature of community and population-level responses to temperature, and indicate that climate warming will likely dampen the seasonal turnover of community structure that is characteristic of high latitude invertebrate communities.

high latitudes, warming occurs at a faster rate than the global average (Overland et al., 2020), so population and community dynamics require dedicated monitoring to understand and mitigate the consequences of climate change.
Climate warming may elicit shifts in the rates of spring emergence and dispersal of arthropods as a result of direct effects on the physiology of organisms or indirect effects through earlier snowmelt (Høye & Forchhammer, 2008;Leingärtner et al., 2014). The timing and order in which species appear during the early community assembly phase in spring will in turn affect their competitive ability and abundance throughout their active season (Bokhorst et al., 2016;Clements et al., 2013;Louette & Meester, 2007). In High-Arctic Greenland, for example, abundances of some spider species have declined over the past 18 years in response to rising temperatures and altered timing of snowmelt (Bowden et al., 2018). Since spiders often play an important role as the first predators to emerge or colonise, shifts in their presence, abundance or behaviour are likely to alter community structure during the assembly phase (Dahl et al., 2018). Changes driven by climate warming in the early phases of community assembly might thus dictate the dominance hierarchy of species and determine community composition throughout the active season (Welch & Harwood, 2014).
As global temperatures rise, species with higher thermal optima are favoured over cold-adapted species. There is debate over whether this turnover will manifest in the long term as a reduction in species richness (if species with low thermal optima become locally extinct), an increase in richness (if more species migrate to higher latitudes from warmer climes; Hickling et al., 2006) or maintenance of richness (if cold-adapted species are simply replaced locally by warm-adapted ones; see Hillebrand et al., 2018). Previous research on naturally heated soil systems provides support for the latter, with increasing temperature leading to a decrease in the evenness of epigeal invertebrates, but no change in species richness (Robinson et al., 2018).
However, inter-annual variation in community structure may be driven by dominant species rather than diversity. For example, temporal stability in a shortgrass steppe plant community was positively correlated with the relative abundance of dominant species, while diversity was not a strong regulator of stability (Sasaki & Lauenroth, 2011). Similarly, in boreal stream systems, long-term compositional changes are predicted to be high even without changes in richness, suggesting that species replacement might be the main mechanism causing climateinduced changes in macroinvertebrate assemblages (Mustonen et al., 2018). Other studies on stream invertebrate communities support this finding, indicating that species replacement and altered dominance patterns, rather than changes in species richness, underlie warming-induced changes in seasonal community turnover (Jourdan et al., 2017). However, it is uncertain whether seasonal species replacements and turnover in community structure are equally resistant to the potential changes in species richness due to warming.
Geothermal systems have been identified as natural laboratories for studying the impacts of environmental temperature on naturally occurring aquatic and terrestrial communities (O'Gorman et al., 2014;Sigurdsson et al., 2016;Woodward et al., 2010). The thermal gradients found in these systems help elucidate how multispecies communities respond to temperature change in a heterogeneous landscape, while capturing the complexity of natural ecosystems. Here, we examined how the assembly and early seasonal turnover of invertebrate communities varied along a soil temperature gradient in a geothermal valley in Iceland. Our goal was to determine the extent to which seasonal community turnover is governed by soil temperature, indicating how climate warming will affect seasonal community dynamics in high latitude areas. We analysed spatiotemporal dynamics in community composition as changes in diversity indices (richness, evenness, diversity) and biomass indices (total abundance, mean body mass, total biomass). We hypothesise that (H1) the diversity and biomass of invertebrates will increase with temperature early, though not later in the season (i.e. due to a more rapid progression through the community assembly phase in warmer environments). Consequently, we predict that (H2) the seasonal variability in community structure will decrease with temperature (i.e. due to the longer duration of established communities in warmer environments). If verified, these hypotheses would indicate decreasing seasonal turnover in high latitude invertebrate communities under climate warming.

| Study site and soil temperature
To investigate the effect of temperature on above-ground terrestrial communities, plots bordering 10 streams were selected in the Hengill geothermal valley, Iceland (64.03°N, 21.18°W). Three plots, each separated by approximately 15 m, were randomly designated on the left and right bank of each stream for a total of 60 plots ( Figure 1a). Soil temperature was monitored for each plot over the sampling period (13 May 2015-3 July 2015) using Maxim Integrated DS1921G Thermocron iButton temperature loggers. The loggers collected hourly measurements, both day and night, and were buried 5 cm below the soil surface (with one additional logger cable-tied to a wooden stake at 5 cm above the soil surface to monitor air temperature). Daily mean soil temperatures closely tracked the temporal pattern of daily mean air temperatures throughout the 2-month sampling period (Pearson's r = 0.70-0.83, p < 0.001 for each plot, n = 52 days; Figure 1b). Data collected during a 2018 field survey also showed that soil temperature at 5 cm below-ground is strongly correlated with surface temperature (Pearson's r = 0.96, p < 0.001, n = 40 plots), justifying the use of below-ground soil temperatures as a measure of the temperature experienced by epigeal invertebrates. Due to the differential geothermal heating of the landscape, a natural soil temperature gradient was present across the plots  (Figure 1b), that is, the plots with the coldest temperatures at the beginning of sampling remained the coldest at the end of sampling (rank correlation of daily mean soil temperatures between dates; Spearman's rho = 0.85, p < 0.001). As such, we used mean plot temperatures calculated across the sampling period for all analyses, henceforth referred to simply as plot temperature. Although there are no data on snowmelt at the plot scale in the lead-up to the study, reconnaissance of the field site 1 month prior to sampling indicated that the entire geothermal valley was snow-free, despite widespread snow cover in the surrounding region ( Figure S1). This suggests that warm and cold plots were likely to have experienced snowmelt at a similar time and thus are unlikely to be at very different stages of assembly post-snowmelt.

| Sampling regime
The epigeal invertebrate communities were sampled using five pitfall traps left for 48 hr in each plot, at four time-points during the study: 19th May, 4th June, 23rd June and 5th July (Figure 1a). Note that this sampling period captures springtime emergence through to the early summertime community, rather than the full seasonal dynamics of the arthropod community (which spans approximately April to October). All invertebrates were identified to species level where possible, but often to higher taxonomic groups where species-level identification was unfeasible (e.g. Diptera and Hymenoptera were mostly identified to family level). We will refer to all taxa as species from henceforth. The activity density of invertebrate species in each plot (henceforth 'abundance') was estimated as the number of individuals of a species found in the five pitfall traps; total invertebrate abundance per plot was calculated as the sum of all species abundances. Individual body masses were estimated by measuring the longest linear dimension and converting to dry weight (mg) using length-weight relationships (see table S2 in Robinson et al., 2018).
The mean body mass of each species and the abundance-weighted arithmetic mean body mass (henceforth 'body mass') of the invertebrate community were estimated for each plot. The biomass of each species per plot was estimated as their abundance multiplied by their mean body mass; total biomass per plot was the sum of all species' biomasses. Alpha diversity of the invertebrate community at each plot was assessed as Shannon-Weiner diversity (H′, using the diversity function in the vegan package in r 4.0.2), species richness (S) and Pielou's evenness (J′ = H′/ln [S]). Seasonal change in these community metrics during the active season was defined as temporal variance, and calculated using the var function in R.

| Univariate statistical analyses
Main and interactive effects of temperature (a continuous fixed effect) and sampling time-point (a categorical fixed effect with four levels) on various community metrics of invertebrate diversity and biomass were analysed using univariate linear mixed effect models (LMEM: lme function in the nlme package) in r. All models included plot nested within bank as a random effect to account for potential pseudoreplication (i.e. three plots on the same stream bank). The main effect of temperature on the temporal variance of community metrics was analysed using LMEM with bank as a random effect.

F I G U R E 1
The distribution and temperature of sampling locations in the Hengill geothermal valley. (a) Location of the 20 stream banks sampled in the study ('=' indicates a greater distance than shown on the map) with the five pitfall traps established in each of three plots at every bank shown in the inset (adapted from Robinson et al., 2018). (b) Daily mean soil (blue-red lines) and air temperature (dashed black line) from 13 May to 3 July 2015. Numbers in the key refer to stream codes used in previous publications on the Hengill system (e.g. Woodward et al., 2010), while 'L' and 'R' refer to the left and right bank of each stream respectively Assumptions were checked from the residuals for all models, and a transformation of log(x) or log(1 − x) was used where necessary to meet the assumptions of normality and homogeneity of variance. To ensure our conclusions were robust to the different levels of taxonomic identification employed for different invertebrate groups, all analyses were repeated after considering taxa at the family level.

| Multivariate statistical analyses
We analysed the same data with Hierarchical Modelling of Species Communities (HMSCs; Ovaskainen et al., 2017), a multivariate approach belonging to the class of joint species distribution models (JSDM; Warton et al., 2015). In these analyses, we considered each visit to each plot of each stream bank as a sampling unit, so that there were 240 sampling units in total (i.e. four visits to three plots at 20 stream banks). We included data on species abundance as the response variable in the HMSC analyses (matrix Y of sampling units × species, including only those species which were present in at least five sampling units). Due to the high proportion of absences, we applied a hurdle model, in which we modelled presence-absence and abundance conditional on presence with two separate models. We modelled presence-absence with a probit regression, and abundance conditional on presence with a linear model, where abundance was log-transformed and scaled for each species to zero mean and unit variance. We included soil temperature (a continuous fixed effect) and sampling time-point (a categorical fixed effect with four levels) as environmental covariates (matrix X of sampled locations × covariates). Note that we examined three different covariate structures: Model 0, containing no environmental covariates (i.e. an intercept-only model); Model 1, containing only the main effects of the covariates; and Model 2, containing the main effects and their interaction. We included data on logtransformed mean body mass as species traits (matrix T of species × traits). Finally, we included random effects to account for the structure of the study design (plot nested within bank) and to evaluate residual species-to-species associations at the level of individual visits (sampling time-point).
We fitted all models with the R package Hmsc (Tikhonov et al., 2020). We sampled the posterior distribution with four independent Markov Chain Monte Carlo (MCMC) chains of 375,000 iterations each, out of which we ignored a transient of 125,000 iterations, and thinned the chain by 1,000 to yield 250 posterior samples per chain and hence 1,000 posterior samples in total. We examined if there was support for the effect of temperature on invertebrate occurrences changing during the study period. To do so, we compared Model 0 (which assumes no effect of temperature) with Model 1 (which assumes a uniform effect of temperature) and Model 2 (which assumes a changing effect of temperature through time) in terms of their fit to the data. We compared the models with the Widely Applicable Information Criterion (WAIC) and with twofold crossvalidation over the 240 sampling units in terms of area under the curve (AUC). We followed Ovaskainen et al., (2017)

and Ovaskainen
and Abrego (2020) by partitioning the explained variation for each species among the fixed and random effects included in the model.

| RE SULTS
For our univariate analyses, the interaction between temperature and time-point was statistically significant for all community metrics apart from body mass ( Table 1). The significant interactive effect of temperature and time-point on diversity was explained by warmer plots having higher diversity than colder plots at the first sampling point, but lower diversity at the last sampling point ( Figure 2a). This pattern was also reflected in the negative effect of temperature on the temporal variance in diversity (Table 1; Figure 3a), meaning the change in diversity through time was greater in cold plots than in warm ones. We observed a significant interactive effect of temperature and time-point on species richness (Table 1; Figure 2b), which was manifested as a significant negative effect of temperature on temporal variance (Table 1; Figure 3b), thus indicating that warmer plots were more stable in terms of richness throughout the active season. The significant interactive effect of temperature and time-point on evenness was explained by warmer plots having higher evenness than colder plots at the first three sampling time-points, but lower evenness at TA B L E 1 Main and interactive effects of temperature and sampling time-point on invertebrate community metrics, and main effect of temperature on the temporal variance of community metrics, tested using linear mixed effects models. Note that 'body mass' is the average size of a species in the community, whereas 'biomass' is the total biomass of all species in the community (body mass × abundance). F and p indicate F-statistics and p-values respectively; p < 0.05 shown in bold  Figure 3c).
There was a significant interactive effect of temperature and time-point on total invertebrate abundance, as warmer plots had higher abundance than colder plots at the last time-point only (Table 1; Figure 2d), but temporal variance was not significantly as- was as a consequence of warmer plots having higher biomass than colder plots at the first and last sampling points, but not at the two sampling points in between (Table 1; Figure 2f). This pattern was also manifested as a significant increase in the temporal variance of biomass with increasing temperature (Table 1; Figure 3f). Similar patterns were detected for all univariate analyses considering taxa at the family level (Table S1; Figures S3 and S4).
For the HMSC presence-absence models, WAIC and predictive power values clearly supported the models that included temperature and time-point (i.e. Models 1 and 2) over the intercept-only model (i.e. Model 0; see Table 2). Including the temperature × time- Each data point represents a single measurement per plot, per sampling timepoint. Solid lines illustrate the predicted fit from significant linear mixed effect models (see Table 1 for full statistical outputs) a negative response and only 8% exhibited a positive response ( Figure S6). The response of species to temperature depended on their body mass, with larger species dominating cool locations and smaller species dominating warm locations (Figure 4a). Body mass was also related to sampling time-point, with smaller species dominating early in the season (Figure 4b). The HMSC abundance model followed a similar trend, with the model including temperature and time-point but not their interaction (i.e. Model 1) as the most parsimonious, albeit with much weaker predictive power than the presence-absence model (r 2 = 0.08; Table S2; Figures S7 and S8).

F I G U R E 3
Effects of temperature on the seasonal change (estimated as temporal variance) of invertebrate community metrics during the active season: (a) diversity; (b) richness; (c) evenness; (d) abundance; (e) body mass; and (f) biomass. Each data point represents the variance of the four sampling time-points. Solid lines illustrate the predicted fit from significant linear mixed effect models (see Table 1 for full statistical outputs)

TA B L E 2 Comparison of the presence-absence HMSC models.
The models either exclude sampling time-point and temperature (Model 0), include their main effects (Model 1) or include their main and interactive effects (Model 2). Explanatory and predictive powers are measured by AUC, where explanatory power is based on comparing the data to predictions of a model fitted to all data and predictive power is based on twofold cross-validation  Koltz et al., 2018;Koltz & Wright, 2020) are now needed to test the generality of these findings for high latitude ecosystems.
During the initial community assembly phase (i.e. the first sampling time-point in mid-May), invertebrate diversity (including both richness and evenness) was higher in warmer plots compared with colder ones, offering partial support for H1 (Figure 2a-c). These findings likely arose from higher soil temperatures in spring promoting earlier snowmelt and enabling invertebrates to colonise or emerge from hibernation or diapause earlier in those plots. For example, an earlier snowmelt treatment in an alpine study led to earlier emergence of arthropods, and an increase in the abundance of predators such as Aranea and Hymenoptera (Leingärtner et al., 2014).
There was no effect of temperature on total abundance at the first sampling time-point here, however, indicating that warming instead may structure communities through changes in dominance hierarchies and relative abundances of species, as observed in High-Arctic plants (Kapfer & Grytnes, 2017), flies (Loboda et al., 2018) and spiders (Bowden et al., 2018).
Earlier snowmelt is also likely to advance plant growth, providing invertebrates in warm plots with more resources (Koltz et al., 2018). Some species' responses to increased availability of resources may be greater than to local environmental conditions (Kankaanpää et al., 2018). For example, the response of the ensign scale insect, Arctorthezia cataphracta, to a combination of nutrient addition and warming was stronger than to warming alone in a 4-year study in Norway (Hågvar & Klanderud, 2009). We found that the random effect of plot explained the vast majority of the variance in the presence of A. cataphracta ( Figure S5), suggesting that an unknown plot characteristic, possibly resource related, may be more important in determining its presence than temperature. Although we are unable to judge the extent to which changes in the composition of our invertebrate community were driven by temperature-mediated effects on vegetation, altered temperature and snow cover regimes can generate transient differences in plant community assembly (Høye et al., 2007;Wipf et al., 2009), leaving legacy effects on the associated invertebrate community.
For instance, plant-plant and plant-invertebrate interactions were very strong drivers of invertebrate community changes in response to alterations in snow cover regimes in an alpine study (Lord et al., 2018). These effects may balance out after the initial community assembly phase as the season progresses, however, with a study in Svalbard indicating that the positive effect of increased vegetation cover on spider abundance diminished at later dates (Dahl et al., 2018).
The species inhabiting cold plots in early spring are likely to be those in the regional species pool with lower thermal optima, but could also be those better able to disperse throughout the landscape. For example, the high mobility and plasticity of life-history strategies in the wolf spider Pardosa palustris make it less susceptible to changes in temperature regimes compared with other arthropod species (Beckers et al., 2020;Hein et al., 2018Hein et al., , 2019. As such, it might colonise cold patches early in spring to avoid competition by other predators (e.g. Pirata piraticus) in warm plots (Robinson et al., 2018). Indeed, we found that P. piraticus responded positively to warming, whereas P. palustris responded negatively ( Figures S6   and S8), indicating that they were colonising disparate habitat patches. Mobile species could therefore be at an advantage in future warmer climates, as they are able to move to potential nano-refugia in heterogeneous environments, but cold-tolerant species might still be lost from the regional species pool if there is not enough habitat heterogeneity to support them (e.g. Alatalo et al., 2017). Soil temperatures in both warm and cold plots in our study closely followed air temperatures throughout the active season, suggesting that cold micro-habitats might disappear in high latitude landscapes under future climate warming.
Abundances were consistently low across the temperature gradient during the first sampling time-point in May ( Figure S2d).
Previous research indicates that cumulative degree-days better predict arthropod abundance than a combination of date and temperature (Bolduc et al., 2013;Schekkerman et al., 2004), offering a possible explanation for why no clear effects of temperature on abundance or biomass were evident during community assembly ( Figure 2). Higher spring temperatures will have a marked effect on species emergence and dispersal (Forrest, 2016), and thus on early season diversity and evenness, but our results indicate that effects on abundance and biomass production will take longer to manifest. Moreover, sampling time-point explained a larger proportion of the variance in the presence-absence of species than temperature ( Figure S5), indicating that cues such as photoperiod may still have an important role to play in moderating the impacts of climate warming on phenology and thus population dynamics (Forrest, 2016).
Examining the seasonal change in diversity showed that warm plots are much more stable across the active period compared with cold plots (though only in terms of Shannon diversity and species richness), partially supporting H2 (Figure 3a,b). Thus, it would appear that warm plots are able to reach their peak productivity earlier than cold plots and maintain a constant species diversity throughout the active season. Further research capturing the entire active season (from the moment of snowmelt to the return of winter snow) is needed to determine whether these dampened seasonal effects at higher soil temperatures were due to faster plant growth and greater dispersal of invertebrates, or earlier snowmelt such that warmer plots were at a more advanced stage of community assembly compared to cold plots. At the peak of the active season (in July), the positive effect of temperature on diversity and evenness that was observed during community assembly became inverted, with diversity and evenness now higher in colder plots (Figure 2a,c). This is supported by our previous findings from Hengill in late summer (August 2012 and July 2013), where both diversity and evenness of invertebrates decreased with increasing temperature (Robinson et al., 2018). This highlights the seasonal context dependency of temperature effects on diversity and the value of temporal replication in sampling regimes to capture the complete story.
It is likely that after the initial community assembly phase, as soil temperature in cold plots increased, more species from the regional species pool were able to colonise these plots. This would suggest that as warm-adapted species are acquired, cold-adapted species are generally not lost from the cold plots. In support of this, we found that many of the taxa in colder plots during community assembly in May (e.g. Entomobryomorpha, Lycosidae, Chironomidae and Linyphiidae) were still present at these plots during the last sampling time-point in July. This indicates that higher temperatures can restructure the community via introduction of new species and changes in relative abundances, rather than through the loss of coldadapted species or species replacement. This may have a significant impact on ecosystem functioning in the light of climate change, as dominant species, rather than species richness, drive ecosystem service delivery (Winfree et al., 2015).
There was greater seasonal change in biomass in the warmer plots, explained by the positive effect of temperature on biomass at the last sampling point (Figure 2f). This is in line with previous research indicating higher body mass of arthropods in sites with warmer summers (Bolduc et al., 2013), and may be due to the presence of larger species at later sampling times ( Figure 4b) and a higher overall abundance in warmer habitats later in the season ( Figure 2d). The general decrease in body mass with temperature ( Figure 4a) does not contradict the higher community-wide biomass at higher temperatures by the end of the season (Figure 2f), which results from a numerically higher abundance of the fewer large species present. The positive effect of temperature on total invertebrate abundance at the last sampling point also follows previous findings from Hengill (Robinson et al., 2018) and elsewhere (e.g. Leingärtner et al., 2014;Loboda et al., 2018). Higher primary production as a consequence of increased cumulative temperature may cause an increase in herbivore biomass (e.g. De Sassi et al., 2012) and consequently also attract more predators at warmer sites.
Although we did not examine patterns in primary productivity here, we found some support for this hypothesis through an increase in the number of Collembola from the largely herbivorous order Poduromorpha, and a greater number of predatory carabid beetles (e.g. Pterostichus diligens and Pterostichus nigrita) and wolf spiders (e.g. Pirata piraticus) at warmer plots ( Figure S8). Moreover, many species which emerge early in the season may extend their period of activity (e.g. aphids; Bell et al., 2015), or increase their voltinism There was a greater prevalence of smaller species as soil temperature increased (Figure 4a), supporting previous findings from late summer in the Hengill terrestrial system (Robinson et al., 2018).
These findings contribute to widespread empirical support for declines in body size with warming across a multitude of taxa and environments (e.g. Sheridan & Bickford, 2011 and references therein, Bowden et al., 2015;Tseng et al., 2018), rooted in temperaturedependent metabolic processes (Brose et al., 2012). Moreover, with an increase in the length of the active season in warmer plots, multivoltine species may be able to maximise their fitness by maturing earlier at a smaller size (Horne et al., 2015). However, body size responses to temperature are highly inconsistent both within and across taxa (see Gardner et al., 2011;Sheridan & Bickford, 2011 for reviews) owing to the heterogeneity in species thermal requirements. Indeed, there was no overall effect of temperature on the abundance-weighted mean body mass of the invertebrate community here (Figure 2e), due to the fact that larger species (such as the predatory P. diligens, P. nigrita and P. piraticus) were more abundant in warmer plots ( Figure S8). The average body mass of the community also increased throughout the season ( Figure S2e), with smaller species found during the first time-point in May than on later sampling occasions (Figure 4b). This is consistent with long-term trends from Siberia, where small arthropods emerged earlier than large ones (Tulp & Schekkerman, 2008), and illustrates the importance of body size for early season structuring of high-latitude terrestrial communities.

| CON CLUS IONS
From our results, we conclude that climate warming will likely dampen the seasonal turnover of community composition that currently characterises high latitude invertebrate communities.
However, while increasing temperature dampens seasonal fluctuations in community diversity, it seems to amplify the temporal variability of community biomass. These effects were driven by contrasting temperature-diversity and temperature-biomass relationships over the course of the summer season. This highlights the importance of seasonal context for a more complete understanding of temperature effects on community structure in climate change research. Microclimate conditions vary along the soil surface and within vegetation, especially in polar terrestrial ecosystems (Convey et al., 2018), helping to maintain suitable conditions for species with a wide range of thermal tolerances in the regional species pool. In this context, our results on the effects of spatiotemporal variation in temperature on an invertebrate community have wider significance for understanding the effects of warming on the structure and dynamics of cold climate invertebrate communities.

ACK N OWLED G EM ENTS
The

DATA AVA I L A B I L I T Y S TAT E M E N T
The data supporting this article are archived with the NERC Environ-