Contrasting latitudinal patterns in diversity and stability in a high-latitude species-rich moth community

Aim: Biodiversity is currently undergoing rapid restructuring across the globe. However, the nature of biodiversity change is not well understood, as community-level changes may hide differential responses in individual population trajectories. Here, we quantify spatio-temporal community and stability dynamics using a long-term high-quality moth monitoring dataset. Location


| INTRODUC TI ON
Determining how ecological communities are structured and how they will respond to global environmental change remains a fundamental challenge for ecologists and biodiversity science (McGill, Dornelas, Gotelli, & Magurran, 2015;Magurran et al., 2018).The Anthropocene is characterized by multiple drivers of change that can impact ecosystems simultaneously, potentially causing heightened pressures on biodiversity and the services that it provides (IPBES, 2019;IPCC, 2014;McGill et al., 2015).Such pressures can prompt rapid and widespread biodiversity responses, which can ultimately lead to destruction of ecosystems and loss of biodiversity (IPBES, 2019).Importantly, current drivers of biodiversity change can impact different populations and taxa differently, resulting in highly variable responses at the community level.For instance, given that thermal niche breadths and sensitivities vary between species, climate change can either reduce or impose constraints on species fitness, abundance and distributions (Bates et al., 2014;Deutsch et al., 2008;Pinsky, Eikeset, McCauley, Payne, & Sunday, 2019;Sunday, Bates, & Dulvy, 2011).Therefore, climate change is driving an accelerated reorganization of ecological communities in space and time (Burrows et al., 2011;Parmesan, 2006;Parmesan & Yohe, 2003;Rosenzweig et al., 2008;Rosenzweig & Neofotis, 2013).
On the other hand, some degree of warming can be expected to promote abundance increases for resident species, given that organisms generally exhibit broader thermal tolerances with increasing latitude (Deutsch et al., 2008;Sunday et al., 2011).Thus, resident species may experience fitness enhancement as thermal conditions get closer to their physiological optima (Bates et al., 2014;Brown, Gillooly, Allen, Savage, & West, 2004).However, how community structure might be responding to these complex changes in species distribution and abundance over space and time remains a critical knowledge gap.This applies in particular to high-latitude systems, given the pervasive paucity of biodiversity monitoring data from these regions (Gillespie et al., 2020;Meyer, Kreft, Guralnick, & Jetz, 2015).
Given divergent responses among different taxa and populations, and across regions, different biodiversity dimensions might exhibit different responses over time.For instance, changes in species richness do not necessarily reflect changes in total abundance or species composition (Dornelas et al., 2014;Hillebrand et al., 2017;Magurran et al., 2018;Schipper et al., 2016).How such discordant changes affect community temporal variability remains poorly established.For many decades (MacArthur, 1955;May, 1973;Tilman, 1996;Tilman, Isbell, & Cowles, 2014), the mechanisms maintaining community stability and its relationship with biodiversity have remained key foci of ecological research.Interest in this area has been accentuated as ecosystems face increasingly interacting anthropogenic pressures (Blüthgen et al., 2016;Craven et al., 2018;Donohue et al., 2016;IPBES, 2019;Isbell et al., 2015;Ma et al., 2017).
Each of these studies to date has focused on a particular level of organization.Of key importance then is the realization that biodiversity change can be assessed at different hierarchical levels: communities consist of species, each of which occurs as populations (then composed of individuals and genes).We may therefore describe community structure and temporal variability, and how those change at different levels, by quantifying change at each level: in species richness, in community composition and/or population dynamics.Each adds a separate dimension, enabling a more comprehensive understanding of how biodiversity is changing over space and time (Hillebrand et al., 2017;Magurran et al., 2018).Yet, we still lack a comprehensive assessment of the combined effects of changes in species diversity and abundance (and further of changes in evenness or dominance; e.g., Hillebrand, Bennett, & Cadotte, 2008;Jones, Collins, Blair, Smith, & Knapp, 2016;Maestre, Castillo-Monroy, Bowker, & Ochoa-Hueso, 2012;Thibaut & Connolly, 2013) and of their effects on community variability across space.For instance, species synchrony, rather than diversity alone, has been shown to strongly influence community stability in plant and animal communities (Blüthgen et al., 2016).Such assessment is urgently needed for insect communities, where recent reports point to widespread (Biesmeijer et al., 2006;Lister & Garcia, 2018;Potts et al., 2010;Sánchez-Bayo & Wyckhuys, 2019;Seibold et al., 2019) but inadequately documented changes (Leather, 2018;Montgomery et al., 2020;Thomas, Jones, & Hartley, 2019).
Here, we quantify the spatio-temporal dynamics of change in total abundance, species richness and community temporal variability in a long-term dataset of systematic moth monitoring (Leinonen et al., 2016(Leinonen et al., , 2017)).Moth communities were sampled continuously between 1993 and 2012 across Finland, spanning c. 10° of latitude (c.1,100 km) from a hemi-boreal to a north-boreal (subarctic) bioclimatic zone.We provide a comprehensive evaluation of temporal changes in community structure and stability, and assess whether these temporal changes are occurring differently along latitudinal or longitudinal gradients.Our findings improve our understanding of how different biodiversity dimensions are changing over space and time, with strong implications for biodiversity change in a species-rich component of high-latitude ecosystems.

| Moth monitoring data
The National Moth Monitoring scheme is coordinated by the Finnish Environment Institute (SYKE), running annually since 1993 (Leinonen et al., 2016(Leinonen et al., , 2017)).Nocturnal moths were sampled using "Jalas" light traps (Jalas, 1960) equipped with Hg vapour bulbs, located mainly in forested areas across Finland (66% of the traps).The proportion of traps located in forest versus other habitats was similar across the bioclimatic zones (Supporting Information Figure S1).The traps were located in the same location from year to year and usually emptied weekly.Sampling occurred every night from early spring (after snowmelt) to late autumn (night frosts and at least one empty trapping period, i.e., 1 week).Sampling effort was constant across years for each trap, but given that sampling aimed to cover the entire moth activity period at each location, the trapping period was naturally longer in more southern traps.Voluntary observers identified the specimens (Leinonen et al., 2016), with a variable number of traps being sampled per year (average of 50 traps per year; Supporting Information Figure S2).The taxonomic skills of the volunteer lepidopterists were typically excellent, and data quality control and cross-checking was carried out by the monitoring coordinators (Pöyry et al., 2011).
The data used for the present study consisted of records of abundance per species collected from 65 traps (unique Trap Code, with unique coordinates) with ≥ 8 years of sampling between 1993 and 2012, yielding 1,000 trap × year combinations.There was no distinct pattern in the number of years sampled in each trap versus latitude (slope = −0.138,t = −0.602,p = .549).We used data sampled between April and October (sparse records outside this period were excluded).For 52 trap × year combinations, traps had not been deployed for the full season, but the results remained unchanged when this minority of data points were removed from the analyses (Supporting Information Supplementary Results I).The data included observations of more than four million individuals and 782 taxa, and specifically, 734 moth species.For the species richness and stability analysis, we excluded records not identified to the species level, but such records were included for measures of abundance.
For each trap × year combination, we pooled the data across the season and calculated abundance per species, total abundance and species richness.We also calculated dominance, as the number of individuals of the most abundant species regardless of species identity, to elucidate further the mechanisms underlying changes in community structure.
Theoretically, dominance patterns could be affected by changes in the identity of the dominant species.If initially dominant species were replaced by incoming species of initially lower abundances, then a potential lag between colonization and subsequent peak abundance might be reflected in a decline in dominance.To test for such imprints of a potential lag phase, for each trap × year combination we regressed the maximum abundance values against the time the dominant species had been present in that particular trap (Supporting Information Supplementary Results II).We fitted a linear mixed-effects model using lmer() from the R package lme4 (Bates, Maechler, Bolker, & Walker, 2015).To account for the fact that some species tend to be inherently more abundant than others over time (i.e., have a higher "carrying capacity"), we included species as a random effect.Overall, the dominant species were species that had been present from the beginning of trapping at the site rather than newly colonizing species.

| Biodiversity metrics and temporal change
For the diversity metrics, we modelled changes in the logarithm of richness and the logarithm of abundance as a function of year, latitude and longitude, and included the three-way interaction between the main effects.We used log 10 -transformed responses to account for pronounced differences in species richness and number of individuals along the latitudinal gradient.All explanatory variables were scaled to a mean of zero and a standard deviation of one.We used lmer() from the R package lme4 (Bates et al., 2015) to fit linear mixed-effects models (LMMs), having trap as random effect to account for repeated sampling, and allowing slopes for year in the random structure to vary for all the models fitted.Specifically, the different models run were: (1) full model with three-way interactions; (2) model excluding interaction terms; and (3) model excluding longitude altogether.We used likelihood ratio tests via the anova() function to compare these models.We used the R package MuMIn (Barton, 2018) to calculate pseudo-R 2 values and the R package ggeffects (Lüdecke, 2018) to plot the marginal effects of the best model for each diversity metric.
We quantified temporal stability in terms of both species richness and aggregated species abundance metrics.First, we calculated the coefficient of variation (CV = σ/μ) of the total number of species found in each trap over time, where smaller values represent greater stability.Note that for a Poisson distribution, μ = σ; therefore, we expected a constant CV = 1 across different values of μ.
Additionally, we calculated community stability as defined by Tilman (1999), that is, the stability of summed species abundances, using the function community_stability() in the R package codyn (v.2.0.0;Hallett et al., 2016).This function aggregates species abundances within replicate and time period and calculates community stability as the temporal mean divided by the temporal standard deviation (μ/σ), where larger values indicate greater temporal stability, that is, lower temporal variation around the mean.We used ordinary least squares regression (OLS) to test for the effect of latitude on community stability metrics (using the logarithm of CV).Because there was heteroscedasticity in the CV response (studentized Breusch-Pagan test: p = .009;lmtest package; Zeileis & Hothorn, 2002), we further assessed the robustness of the OLS estimates via regression with robust standard errors, using the package sandwich (v.2.5-0; Zeileis, 2004).Given that the overall direction and magnitude of the pattern were consistent (only wider confidence intervals for the parameter estimates in the latter), we report OLS results only.We evaluated the potential effects of calculating stability metrics when sampling did not occur over consecutive years by running the same analysis excluding four traps that had temporal gaps > 1 year over the sampling period.Finally, to examine how individual species patterns contributed to overall stability at the community-level, we calculated the index of community synchrony proposed by Loreau and de Mazancourt (2008), using the function synchrony() in codyn.This metric compares the variance of aggregated species abundances with the summed variances of individual species, and can deal with variable species richness values.To test the effect of latitude, we fitted a linear model to the logit transformation of synchrony values, and further assessed the relationship between stability and synchrony with OLS as well.All analyses were conducted in R (R Core Team, 2017).

| RE SULTS
Overall, species richness increased during the monitoring period.
This effect was more pronounced with increasing latitude (steeper positive rates), while longitude had no effect on species richness change (Figures 1a and 2a; Table 1; Supporting Information Figure S3).The most parsimonious model included year, latitude and their interaction (Table 1; Supporting Information Table S1; longitude effects were not significant [ 2 models 1−2(d.f.=3) = 0.224, p = .974; 2 models 2−3(d.f.=1) = 3.682, p = .055]).We further detected an overall decline in total abundance (no spatial interactions; Figures 1b and   2b; Table 1; Supporting Information Table S2; [ 2 models 1−2(d.f.=4) = 3.092, p = .543; 2 models 2−3(d.f.=1) = 0.378, p = .539]),indicating that despite the increase in richness, the number of individuals in the community decreased over the monitoring period.The decline in total abundance was also evident in the numerically dominant species (Figure 1c).Furthermore, we showed that the dominant species in each trap and year was mostly present throughout the period sampled in each trap (Supporting Information Supplementary Results II).
On average, the moth community experienced increased temporal variability with increasing latitude, and this applied to both species richness and aggregated population trends metrics.CV increased with latitude, indicating that the variability in relation to the mean species richness increased poleward (Figure 3a; Table 2).
A similar pattern emerged in terms of the summed species abundances stability metric, where stability likewise decreased with increasing latitude (Figure 3b; Table 2).These results remained unchanged when excluding traps that had temporal gaps > 1 year over the sampling period (Supporting Information Table S3).The synchrony in population trends among species increased strongly poleward (Figure 4a; Table 2).Finally, there was a strong negative relationship between community stability and synchrony (Figure 4b), indicating that higher asynchrony was associated with increased stability, which occurred predominantly across traps at low to midlatitudes within our data.

| D ISCUSS I ON
While the impacts of global change on specific levels of community organization are relatively well characterized, their relationship to changes across different organizational levels from populations to communities are less so.In this study, we uncovered contrasting temporal patterns of change in the structure and stability of a highlatitude insect community.Against a backdrop of increasing species richness over time, we found a clear latitudinal pattern where community dynamics tended to be less stable the higher the latitude.
These patterns were underpinned by faster increases in species richness towards the pole, while total abundance consistently declined, mainly associated with a decline in the most dominant species.Our findings reveal important temporal latitudinal patterns in this boreal community.Below, we discuss each finding in turn.

| Spatio-temporal changes in community structure
The overall increase in species richness reported here in Finnish moth communities is consistent with widespread reports of poleward shifts in species' geographical ranges, including terrestrial invertebrates and Lepidoptera (Chen et al., 2011;Parmesan & Yohe, 2003;Pöyry et al., 2009).It is also in line with previous reports suggesting that the species composition of Finnish moths has, on average, become "south-westernized", with increasing prevalence of species historically found at lower latitudes (Leinonen et al., 2016).
What we add to this picture is that this effect has unfolded proportionally faster poleward-a pattern consistent with asymmetric latitudinal range shifts, where the leading edge advances faster than the trailing edge (Sunday, Bates, & Dulvy, 2012).A matching pattern has been reported for Finnish birds, where nothern species have shifted their ranges poleward faster than southern species (Virkkala & Lehikoinen, 2014).
In our study, the asymmetrical increase in species richness is likely underpinned by faster rates of warming occurring at high latitudes (Loarie et al., 2009;IPCC, 2014), as well as within Finland (Virkkala & Lehikoinen, 2014).Therefore, populations at the poleward range margin may be able to expand quickly into new geographical areas, following faster shifting local conditions (Schmid, Dallo, & Guillaume, 2019;Sunday et al., 2012).On the other hand, trailing-edge populations at the equatorial range margin, which are likely expericencing slower local changes, may not be moving as quickly, nor may they be facing conditions beyond their thermal optima.Importantly, most Finnish moth species have a southern edge margin extending far south of our study region and will thus be far from their upper thermal limits (for distribution maps of all European moth species, see e.g., www.lepid optera.eu).Furthemore, how such range shifts unfold likely depends on species traits, such as dispersal ability, and on the availability of suitable habitats (Pöyry et al., 2009;Warren et al., 2001).

| Contrasting patterns in different metrics
In rapidly warming high-latitude regions, the number of individuals and species is expected to increase (e.g., Deutsch et al., 2008;Elmhagen et al., 2015).While we found strong support for the latter, we found the opposite for total abundance.Matching recent reports of significant declines among insects (Biesmeijer et al., 2006;Dirzo et al., 2014;Montgomery et al., 2020;Potts et al., 2010;Sánchez-Bayo & Wyckhuys, 2019;Seibold et al., 2019), we showed a substantial decline in the overall abundance of moths, which was associated with declines of the most abundant species, using long-term systematic monitoring data.Given the fundamental role of insects in driving ecosystems structure and function, such abundance declines are likely to negatively impact
Our results highlight that changes in richness and abundance can be decoupled from each other (Schipper et al., 2016).In a similar vein, Antão et al. (2019) reported how temperature-related trends in the number of individuals and number of species were not only decoupled, but exhibited geographically contrasting patterns during several decades of climate change.Such contrasting trends between richness and abundance might be explained by the ability of species to exploit microrefugia (Suggitt et al., 2018), allowing populations to persist in smaller pockets of suitable climate in the face of environmental changes (Aalto, Riihimäki, Meineri, Hylander, & Luoto, 2017;Greiser, Meineri, Luoto, Ehrlén, & Hylander, 2018;Maclean, Suggitt, Wilson, Duffy, & Bennie, 2017;Meineri & Hylander, 2017).
A retraction to smaller parts of the landscape is likely to reflect into a reduction in population size, without necessarily causing a matching effect on species occurrence.
Incidentally, we note that the patterns observed are the opposite to those expected from observation biases, where a decreased abundance should result in fewer species crossing a potential observation threshold, and thus in reduced (not increased) species richness.Thus, our results offer strong evidence that a one-eyed focus on a metric such as species richness (as reflecting occupancy) could be masking widespread declines in abundance (reflecting population-level trajectories).Such abundance declines could be linked to several underlying factors, such as changes in temperature and precipitation patterns resulting from ongoing climate change (IPCC, 2014) and phenological mismatches with the availability of food resources (Renner & Zohner, 2018;Schmidt et al., 2016;Visser & Both, 2005;Visser & Holleman, 2001).They are also likely to be coupled to habitat degradation or reduction (Fox et al., 2014;Habel et al., 2019;Potts et al., 2010).Evaluating different biodiversity metrics reveals different and complementary processes affecting communities over time.Our findings reinforce the notion that changes in one biodiversity dimension fail to represent the suite of restructuring trajectories occurring over space and time.

| Shifting patterns of dominance
The overall decline in abundance was associated with a decline of the quantitatively most dominant species (the top 10 most Responses to the combined effects of year and latitude of: (a) species richness; and (b) total abundance, showing marginal effects corresponding to the best models fitted to each metric according to Table 1.Log S and Log N refer to the logarithm of species richness and of number of individuals, respectively.
For each bioclimatic zone, we show the response at the level of the mean latitude value (colours as in Figure 1).Note that year and latitude were scaled to a mean of zero and unit variance [Colour figure can be viewed at wileyonlinelibrary.com]Note: For the model selection results, see the Supporting Information (Tables S1 and S2, respectively).The estimated coefficients are shown with the 95% confidence intervals, and marginal (fixed) and conditional (both fixed and random effects) pseudo-R 2 values.everywhere (e.g., Inger et al., 2015;Jones & Magurran, 2018).Unlike the pattern of the dominant species accounting for roughly the same proportion of the assemblage abundance revealed by a recent metaanalysis across taxa and ecosystems (Jones & Magurran, 2018), we found an overall decline in dominance.We further showed that this decline is not underpinned by lags between colonization and peak abundance.Overall, our findings are consistent with faster colonizations and influx of species that are likely initially rare (prompting the increase in species richness), while abundance declines are mostly driven by the resident dominant species, which are mostly present throughout the sampling period.However, given that we did not track species identity, we are unable to further identify the potential effects of changing dominant species identity and traits, which fundamentally determine ecosystem functioning and stability (Grman, Lau, Schoolmaster, & Gross, 2010;Hillebrand et al., 2008).
The combined effects of temperature change (likely prompting the increase in richness), coupled with other drivers of change, might be impacting the abundance of populations in a stronger adverse way.
One such driver shown to affect moth communities in Finland is soil nutrient enrichment caused by nitrogen deposition (Pöyry et al., 2017).Furthermore, Leinonen et al. (2016) noted that the community evenness (as quantified by Fisher's α index) has increased between 1993 and 2012, suggesting shifts in the species abundance distribution consistent with declines in the most abundant species and/ or increases in species richness.Declines in the dominant species can potentially lead to disruption of community structure, species interactions and ecosystem services, in particular because common species may contribute disproportionately to such processes (Gaston, 2010;Hillebrand et al., 2008;Winfree, Fox, Williams, Reilly, & Cariveau, 2015).In the case of insects, such effects are likely to propagate to higher trophic levels.

| A latitudinal pattern in temporal variability
For half a century, the question of how much communities with different characteristics vary over time has been one of the most researched topics of ecology (Craven et al., 2018;Isbell et al., 2015;MacArthur, 1955;May, 1973;Tilman, 1996;Tilman et al., 2014).The vast majority of studies on the effects of biodiversity on stability have largely focused on loss of species (Tilman et al., 2014), whereas species asynchrony may be a stronger driver of stability (Blüthgen et al., 2016).Furthermore, most research on stability has focused on plant communities (Donohue et al., 2016).In insects, models focusing on dynamic equilibrium between disturbances and stability might be more suitable for predicting diversity patterns than models focusing on interspecific competition and resource partitioning (e.g., Levin & Paine, 1974).For instance, plant species richness is strongly affected by disturbance dynamics, and given that moths are herbivores, their richness is strongly correlated with plant richness (Haddad et al., 2009;Schuldt et al., 2019;Siemann, Tilman, Haarstad, & Ritchie, 1998).In this study, we relate the effects of an overall increase in moth species richness across a wide latitudinal range to regional patterns in community variability.We show that such an increase was associated with higher variability with increasing latitude.
Overall, variability over time and synchrony clearly increased towards the North.Three other latitudinal patterns co-occur with these trends: species gain is occurring proportionally faster poleward; abundances decrease; and some northern moth populations display recurrent abundance outbreaks (Jepsen et al., 2008;Klemola, Huitu, & Ruohomäki, 2006;Leinonen et al., 2016;Supporting Information Figure S4).Although the same species are present in the South of our study area, abundance outbreaks are rarer and milder there (Klemola et al., 2006).Thus, we are far from being able to causally relate the patterns to any particular driver.What we can do is to point to a similarity with patterns in other organisms.For instance, temporal patterns in population stability have been examined in small mammals, with increasing amplitude and the propensity for cyclicity increasing with latitude (Hansson & Henttonen, 1985, 1988).Elmhagen et al. (2015) report widespread changes in ranges, abundances and cyclic patterns for mammals and birds in boreal forest and alpine tundra.Moreover, the cyclic outbreaks of moth species have been associated with warming-induced range expansions across northern latitudes (Jepsen et al., 2008).Indeed, lepidopteran

F
Temporal patterns of: (a) species richness; (b) total abundance; and (c) dominance (all on a logarithmic scale), plotted for each of the four bioclimatic zones across Finland (from South to North: HB = hemi-boreal; SB = south-boreal; MB = mid-boreal; NB = northboreal).Grey dots represent the number of species or individuals per individual trap in each year, and coloured lines indicate the temporal trend per trap.The y axis was truncated at the lower limit to allow better visualization of patterns across latitude (excluding seven and three traps with very low number of species and individuals in the NB zone, respectively) [Colour figure can be viewed at wileyonlinelibrary.com] Changes in species richness and total abundance as a function of year and latitude Metrics of community temporal stability, shown as a function of latitude: (a) coefficient of variation of species richness (Log CV Sps); and (b) stability.Each dot represents a unique trap, coloured according to bioclimatic zone as in Figure 1; the size of the dots is proportional to the number of years sampled per trap [Colour figure can be viewed at wileyonlinelibrary.com]Results from models of community stability and synchrony as a function of latitude The estimated coefficients are shown with the 95% confidence intervals and R 2 values.*p < .05;**p < .01;***p < .001.F I G U R E 4 (a) Community synchrony as a function of latitude.(b) Community synchrony in relation to community stability.For synchrony, zero indicates perfect asynchrony between species, and one indicates perfect synchrony.Colours indicate the bioclimatic zone as in Figure 1, and the size of each dot is proportional to the number of years sampled per trap [Colour figure can be viewed at wileyonlinelibrary.com] average across each bioclimatic zone are listed in Supporting Information Figure S4).A similar pattern for dominance was described by Yang et al. (2017) for a temperate steppe plant community subjected to warming, where temporal stability was reduced mainly through the decline in abundance of the dominant species.These results suggest the possibility that dominant species are becoming relatively more abundant over time might not hold true communities have recently shown pronounced temporal change, specifically prompted by climate change and with clear latitudinal patterns, including northward range shifts(Pöyry et al., 2009) and increased multivoltinism(Pöyry et al., 2011).Our findings suggest dramatic changes in the structure of a high-latitude moth community across a geographical gradient of 1,100 km, covering climates stretching from near temperate to subarctic.With different patterns emerging in different metrics of community structure, our study highlights the crucial point that increases in species richness may actually mask a substantial loss of individuals.Thereby, our findings suggest that the ongoing species redistribution is systematically affecting community structure, population abundances and stability over time, leading to compounded and partly opposing effects of global change.In turn, these can negatively affect higher trophic levels and ecosystem structure.These contrasting trajectories of different biodiversity dimensions are driving a reshuffling of ecological communities, potentially leading to no-analogue communities (i.e., communities of types not seen before;Williams & Jackson, 2007;Williams, Jackson, & Kutzbach, 2007).With increasing future variation in climatic drivers(IPCC, 2014), potentially more extreme outcomes can be expected, particularly given the combined changes across the different aspects of community structure.A better understanding of the mechanisms underpinning such changes, and their effects on community stability, is paramount in order to preserve biodiversity and the services it provides.Exposing the consequences of global environmental change requires a versatile perspective on community structure to dissect diverse patterns of spatio-temporal dynamics.