Late Quaternary chironomid community structure shaped by rate and magnitude of climate change

Much is known about how climate change impacts ecosystem richness and turnover, but we have less understanding of its influence on ecosystem structures. Here, we use ecological metrics (beta diversity, compositional disorder and network skewness) to quantify the community structural responses of temperature‐sensitive chironomids (Diptera: Chironomidae) during the Late Glacial (14 700–11 700 cal a bp) and Holocene (11 700 cal a bp to present). Analyses demonstrate high turnover (beta diversity) of chironomid composition across both epochs; however, structural metrics stayed relatively intact. Compositional disorder and skewness show greatest structural change in the Younger Dryas, following the rapid, high‐magnitude climate change at the Bølling–Allerød to Younger Dryas transition. There were fewer climate‐related structural changes across the early to mid–late Holocene, where climate change was more gradual and lower in magnitude. The reduced impact on structural metrics could be due to greater functional resilience provided by the wider chironomid community, or to the replacement of same functional‐type taxa in the network structure. These results provide insight into how future rapid climate change may alter chironomid communities and could suggest that while turnover may remain high under a rapidly warming climate, community structural dynamics retain some resilience.


Introduction
Current climate change is occurring at an unprecedented rate that will continue over the coming decades (Smith et al., 2015). The rate and magnitude of climate change affects the capacity of ecosystems to absorb climate shocks (Overpeck et al., 1991;Skelly et al., 2007;Grimm et al., 2013). A fast rate impacts the abilities of organisms to respond effectively and efficiently to changes (Shayegh et al., 2016), while a high variability influences the extent of ecosystem response (Stireman et al., 2005). The nature of community-level response to climate-induced stress will depend on the structural organisation and connectivity of taxa within the ecosystem (Dunne et al., 2002;van Nes and Scheffer, 2005;Scheffer et al., 2012). Structure and connectivity shape both resilience to stress and type of response to the stressor, whether smooth (linear) or abrupt (nonlinear: Burkett et al., 2005;Scheffer et al., 2012). Here we aim to improve our understanding of how community structures respond to different rates and magnitudes of climate change by compositional and network analyses of micro-faunal datasets in lake sediments spanning major climate cycles over the last 15 000 years. The payoff for understanding how community structures have responded to recent climate change is an improved ability to anticipate the conditions for future ecosystem stability.
Climate is thought to drive changes in ecosystems through two mechanisms: directly through physiological stress; and indirectly by changing taxonomic interspecific interactions, and thus the ecosystem structure (Harley, 2011). Climate change is projected to drive taxon richness and biodiversity loss (Thomas et al., 2006;Foden et al., 2013), with the extent of changes in the ecosystem depending on functional diversity (Chapin et al., 2000). Functional diversity means the number and range of functions that taxa have within the ecosystem (Petchey and Gaston, 2006). Thus, the loss of a taxon can have a different effect on ecosystem structure depending on the connectivity and function of the taxon within the ecosystem (Dunne et al., 2002). A number of studies have started to evaluate the effect of losing functional groups of taxa on community structure. For example, Doncaster et al. (2016) and Wang et al. (2019) explored the effect of taxon-type losses from aquatic communities, suggesting that the early loss of less common, weakly interconnected taxa may leave a predominance of more strongly interconnected taxa defining a more rigid structure prone to abrupt collapse. The ability of a taxon to persist within an ecosystem relates to its tolerance levels, with greater magnitudes of stress often driving greater assemblage change (Cao and Hawkins, 2005). Thus, investigating how biodiversity, taxonomic interactions and community structure are influenced by different magnitudes of climate change is fundamental to understanding future ecosystem functioning and stability.
Here, we seek to detect signals of compositional change in response to past large-scale temperature changes. We focus on two periods: the relatively rapid, high-magnitude climate change during the Last Glacial-Interglacial Transition (LGIT), and the more gradual, relatively low-magnitude climate change experienced during the Holocene. The LGIT was a period of rapid, high-amplitude climatic instability (Hoek, 2001). It was composed of two millennial-scale oscillations; the Bølling-Allerød interstadial, a period of abrupt and significant warming (c. 14 700-12 800 cal a BP) (Dansgaard et al., 1993;Buizert et al., 2014;Rosén et al., 2014), and the Younger Dryas stadial, an abrupt return to cool, glacial conditions (c. 12 800-11 700 cal a BP) (Dansgaard et al., 1993;Tarasov and Peltier, 2005). In contrast, the Holocene climate has been relatively stable (Wolff et al., 2010). Here, we divide the Holocene into two sections; the early Holocene (11 700 cal a BP to 8200 cal a BP) and mid-late Holocene (8200 cal a BP to present). The early Holocene (11 700 cal a BP to 8200 cal a BP) experienced initial rapid climatic warming following deglaciation (Birks and Birks, 2008;Aarnes et al., 2011). The mid-late Holocene (8200 cal a BP to present) was less climatically changeable, with weakened orbital forcing during the northern hemisphere summers (Wanner et al., 2008), particularly in high latitude regions (Balascio and Bradley, 2012). While there were a number of widely recognised climatic events during the mid-late Holocene, the timing of such events, e.g. the Holocene Thermal Maximum and neoglaciation, was time transgressive (Kaufman et al., 2004;McKay et al., 2018), variable across high-latitude regions (Renssen et al., 2009(Renssen et al., , 2012Briner et al., 2017;Geirsdóttir et al., 2019), and spanned the mid-late Holocene boundary in some regions (Miller et al., 2005;Salonen et al., 2011;Badding et al., 2013). We therefore group the 'mid' and 'late' Holocene sections together.
We focus on temperature-sensitive chironomid (Diptera: Chironomidae; non-biting midges) records from northern Europe (Norway and Scotland) as indicators of the effect of temperature on aquatic ecosystems. Chironomid records have frequently been used as palaeoecological proxies of past temperature changes (Brooks and Birks, 2001;Brooks, 2006;Axford et al., 2017), thus it is thought that such records may provide a greater insight into the complex relationships between climatic variability and ecosystem response (Willis et al., 2010;Birks et al., 2016). We test for differences in the ecological response to the relatively rapid, highmagnitude climate change at the Bølling-Allerød to Younger Dryas transition and the more gradual, relatively low-magnitude climate change experienced at the early-mid-Holocene transition. We employ a suite of ecological metrics to compare the assemblage composition before and after the climatic change, including ordination, taxon richness, beta diversity, compositional disorder and skewness. These metrics have previously been used to assess chironomid structural change within spatial datasets (Mayfield et al., 2020); here we use them to assess changing chironomid community structure in temporal records for the first time. To assess how the assemblages may have changed under continued Bølling-Allerød and early Holocene climatic conditions, we test for evidence of departures from a status quo with autoregressive integrated moving average (ARIMA) forecasts. We test for a response in chironomid assemblage structure to changing climate by comparing the ecological metrics with the NGRIP 15 000-year isotope record trends and fluctuations. To distinguish true pattern from random noise in ecological datasets, we follow the principles in Telford and Birks (2011) for minimising the likelihood of drawing conclusions from spurious correlations. Our criterion for detecting non-random patterns in empirical datasets is that the observed metric must explain more of the variation in a response than is explained by application of the metric to randomised datasets in at least 95% of 10 000 randomisations. We expect to see a greater response in the ecosystem metrics from the rapid, large magnitude change in the Bølling-Allerød to Younger Dryas if the chironomid assemblages were responding to climate change. For example, if the rapid, abrupt change in climate drives a faster and greater magnitude change in taxon presence, this could have a greater effect on the organisation and functionality of the ecosystem. In contrast, the more gradual, lower-magnitude change between the early and mid-late Holocene should produce smaller-scale changes in the ecosystem metrics. For example, under lower stress conditions, the ecosystem may have greater capacity to adjust with the changing boundary conditions; thus any potential changes in taxonomic composition were likely to have a lesser effect on the ecosystem structure.

Quantifying climate change
Oxygen isotope (δ 18 O‰) records from the NGRIP Greenland ice core were used to provide an independent record of relative air temperature for the recent glacial-interglacial period (Wolff et al., 2010) (Fig. 1, upper panel). To assess the extent of temperature change across the Late Glacial and Holocene epochs, median δ 18 O (‰) values were calculated for the Bølling-Allerød and Younger Dryas, and early and mid-late Holocene periods. Median δ 18 O (‰) values indicate a relatively large difference in temperature between the Bølling-Allerød (-38.61) and Younger Dryas (−40.77). The Holocene, in contrast, has been relatively more climatically stable (Fig. 1, upper panel) (Wolff et al., 2010), with little difference between the median δ 18 O (‰) values for the early (−34.95) and mid-late (−35.07) Holocene subsections. To ascertain whether the magnitude of change between the observed medians exceeded the expectations for random noise, the observed differences in the medians were compared with differences in median metrics obtained from 10 000 randomised replications of the Late Glacial and Holocene isotope record sections. The randomised datasets were created by randomly reordering the isotope record sections. Each randomised dataset was partitioned into the corresponding Bølling-Allerød and Younger Dryas and early and mid-late Holocene subgroups. Median values were calculated for the randomised subsets. Differences between the medians were calculated to produce 10 000 differences in median values. These differences were ranked, along with the observed difference in median values. For the Late Glacial, the observed difference in medians fell in the lower 2.5% of the ranking, indicating that the observed difference had a less than 5% chance of being random occurrence ( Supplementary  Fig. 1). For the Holocene, the observed difference in medians did not fall outside the upper or lower 2.5% of the ranking, indicating that the observed difference for the Holocene isotope medians was not discernible from random occurrence ( Supplementary Fig. 1).
To quantify the different rates and magnitudes of climate change during the Late Glacial and Holocene periods, we calculated the cumulative deviation from the long-term mean (Dugmore et al., 2007). To do this, long-term rolling mean values were calculated over periods of 200 years (10 samples) ( Supplementary Fig. 2). The gradient of the slope indicated the relative velocity of the climate change and the amplitude of the curve indicated the magnitude of the climate change. The greatest changes in the rate and magnitude of the isotope trend occurred at the start of the Bølling-Allerød and Younger Dryas periods, with smaller magnitude fluctuations within each period ( Fig. 1, lower panel). In the Holocene record, the δ 18 O (‰) record showed the greatest difference in cumulative deviation from the long-term mean at the start of the Holocene. The cumulative deviation from the long-term mean gradually decreased through the early Holocene, indicating a slower rate of change. There were rapid, small magnitude changes in the cumulative deviation trend in association with the Preboreal Oscillation, 9.3 ka and 8.2 ka events. There were few substantial cumulative changes in the isotopic trend throughout the mid-late Holocene, with little change in association with the Holocene Thermal Maximum (c. 9-5 ka in Greenland; Kaufman et al., 2004), neoglaciation (c. 4-3 ka in Greenland; Briner et al., 2017), Medieval Climate Anomaly (950-1250 CE; Mann et al., 2009) or Little Ice Age (c. 1400-1700 CE; Mann et al., 2009), indicating relative climatic stability through the mid-late Holocene.

Empirical chironomid datasets
Larval chironomid assemblages were used here as indicators of ecological stability through time from European datasets (Fig. 2). The authors of the chironomid datasets are indicated in Table 1. Taxa were identified to genus or speciesmorphotype level using standardised subfossil taxonomy (Brooks et al. 2007). While there are a large number of chironomid records available spanning Late Glacial and Holocene time periods, the records for this study were selected based on high sampling resolution, consistency of taxonomic resolution, and age-model quality to allow comparison among sites.
The Late Glacial chironomid records used were from Scotland, UK: Ashik and Abernethy (Brooks et al., 2012b) and Muir Park (Brooks et al., 2016). These lakes span a range of environments, including coastal and inland locations, and varying proximities from the Loch Lomond ice limits (Chandler et al., 2019;Lowe et al., 2019). Age models were available for the Ashik and Abernethy records from the original publications; however, no age model was available for the Muir Park record beyond the tephras (with associated dates) found within the stratigraphy. Estimated ages were calculated for the Muir Park record in this paper ( Supplementary Fig. 3), based on the tephrochronology, to establish a broad chronostratigraphy. The age-depth model for the Ashik record was also based solely on tephrochronology, using the Borrobol and Penifiler Figure 1. Oxygen isotope records from the NGRIP record in Greenland spanning the Last Glacial-Interglacial Transition and Holocene (upper panel) (Wolff et al., 2010). Lower isotopic values are indicative of lower temperatures. Dates have been corrected to cal a BP from yrs B2K, following Wolff et al. (2010). Mean isotopic values are indicated by horizontal red lines. The rate of change, calculated as the cumulative deviation from the rolling long-term (200 years) mean (lower panel), indicates a greater rate and magnitude of isotopic change at the onset of the Bølling-Allerød (BA), Younger Dryas (YD) and early Holocene.
The Holocene chironomid records used were from Norway and span a latitudinal gradient, covering arctic, coastal and alpine locations: Horntjernet from Finnmark, Arctic Norway (this paper), Bjornfjelltjønn on the northwest coast of Norway (Brooks, 2006) and Holebudalen in southern Norway (Velle et al., 2005a). Age models for the Bjornfjelltjønn and Holebudalen records were provided by the original authors. The age model for the Horntjernet record was provided by Rijal et al. (2020).
To assess the comparability of the fossil records, detrended correspondence analyses (DCA) were used to plot the records in ecological space, using vegan (Oksanen et al., 2013) in R statistical software v.3.6.0 (R Core Team, 2019). The relative positions of the Bølling-Allerød and Younger Dryas samples and the early and mid-late Holocene samples indicated whether the assemblages were responding similarly to environmental stressors. Timetrack (Simpson and Oksanen, 2019) and ordisurf (Oksanen et al., 2013) were used to passively plot the temporal records onto the ecological space of the Norwegian calibration set (Brooks andBirks, 2001, 2004) to determine the influence of mean July temperature on the fossil chironomid records in R statistical software v.3.6.0 (R Core Team, 2019). To further test whether the Norwegian calibration set was an appropriate fit for the fossil datasets, a goodness of fit analysis was performed using squared residual distance (Juggins and Birks, 2012) (Supplementary Figs. 4 and 5).

Metrics of structural change
Four ecological metrics, taxon richness, beta diversity, compositional disorder and network skewness, were selected as indicators of ecosystem change for the temporal empirical chironomid assemblages. Beta diversity, compositional disorder, and network skewness reflect ecosystem structural change, through taxonomic turnover and taxon organisation (Baselga, 2010;Doncaster et al., 2016;Wang et al., 2019;Mayfield et al., 2020). Palaeoecological datasets often have irregular time series (Glew et al., 2001), including the records analysed here. Our approach in this paper is to sample changes from one climatic period to the next, rather than change over continuous time, so the records were not corrected for temporal resolution.
Taxon richness has long been used as a measure of the number of taxa present in a sample; however, it does not account for taxon abundance, rarity or ecosystem structure (Veech, 2018;Wang et al., 2019). It was used here as a simplistic characterisation of assemblage change over time and during climatic stress. Taxon richness can be affected by specimen preservation and taxonomic resolution (Velle and Larocque, 2008); in an attempt to account for this, the records analysed here were screened for taxonomic consistency. Here, we used rarefaction as a measure of taxon richness to account for variations in count sums when comparing samples. To calculate rarefaction, taxa counts were rounded to integers and the minimum count sum was used to calculate rarefaction for each site, using the rarefy function in the vegan package (Oksanen et al., 2013) in R statistical software v.3.6.0 (R Core Team, 2019).
Beta diversity, a measure of variation in taxonomic composition (or turnover) (Legendre and De Caceres, 2013), was used to indicate the similarity of assemblages in consecutive time slices. Beta diversity was calculated for each sample following the variance partitioning framework of Legendre and De Caceres (2013), using the beta.div function and 'Hellinger' method in the adespatial package (Dray et al., 2019) in R statistical software v.3.6.0 (R Core Team, 2019). The taxon assemblage data were transformed by the 'Hellinger' method within the beta.div function. This method was selected because a diversity value is calculated for each time slice, as opposed to pairwise site comparisons, which are normally used for spatial datasets. Higher beta diversity values indicate greater variation (or turnover) between assemblages. As assemblages gradually change over time, beta diversity should be low and constant, as there is a high degree of taxonomic overlap. Abrupt changes in the ecosystem composition should cause a sudden rise in beta diversity as the assemblage experiences greater taxonomic turnover.
Compositional disorder, a measure of nestedness, can be used as a proxy for unpredictability in community composition. Compositional disorder is expressed as degree of disorder (°disorder) and measured on a scale from 0°(perfectly nested sequences) to 100°(perfectly disordered sequences) using subsets of 10 consecutive samples (Doncaster et al., 2016;Mayfield et al., 2020). A higher value indicates that the sample was less predictable from previous samples, due to greater dissimilarity between the assemblage compositions. Under low stress conditions, taxa should continually arrive and exit the ecosystem in response to fluctuating environmental variables; this should produce a constant level of°disorder with little directional change (Doncaster et al., 2016;Mayfield et al., 2020). Sudden changes in ecosystem composition should produce an abrupt change in°disorder. As taxa rapidly appear, the proportion of common taxa between samples is likely to diminish, causing greater unpredictability between samples and increasing°disorder values. A sudden decline in taxon abundance is more likely to produce a nested assemblage of core, common taxa, producing a decline in°d isorder and signifying greater predictability. Network skewness measures the hierarchical distribution of nodal degree, i.e. the number of interspecific interactions  (Velle et al., 2005a) attributable to each species (termed species degree) (Wang et al., 2019). Network skewness was calculated in MATLAB (ver. R2017b) following the method developed by Wang et al. (2019). This method requires a calibration dataset to identify frequently co-occurring chironomid taxon pairs. In this study, we used the Norwegian calibration dataset (Brooks andBirks, 2001, 2004); goodness of fit analyses confirmed that the Norwegian calibration dataset was an appropriate fit (using squared residual distance (Juggins and Birks, 2012), Supplementary Figs. 4 and 5). The most frequently co-occurring chironomid taxon pairs were identified as the pairs occupying the upper two quartiles of positive values for Cramér's association coefficient (Q2 of V+) in the Norwegian calibration dataset (Brooks andBirks, 2001, 2004). Species degree was calculated for each fossil assemblage based on the co-occurrences of taxa in the Norwegian calibration set. Network skewness was measured as the skewness of the frequency distribution of species degree; see Wang et al. (2019) for further details. A positive skewness value indicates a large proportion of taxa with few interspecific interactions and a small number of taxa with many interactions. Environmental degradation may reduce or change the availability of micro-niches, causing the taxa to reorganise structurally, hence increasing the number of taxon interactions (Albert et al., 2000). This creates a larger proportion of highly connected taxa, producing a less positively skewed distribution of taxonomic interactions (Wang et al., 2019). A prerequisite and limitation of this method is that the taxa in the fossil assemblages must be present in the calibration set. Thus, the chironomid datasets were selected and filtered with this in mind, with some sub-categories of taxa being grouped to genus. These filtered datasets were used for the skewness analyses.
Comparisons between the filtered and full datasets indicated that the filtering process had a limited effect on the taxon richness, beta diversity and compositional disorder analyses (Supplementary Figs. 6 and 7).
Testing structural change across high-and low-magnitude climate change We tested the magnitude of structural change in the chironomid assemblages in the Late Glacial, between the Bølling-Allerød and Younger Dryas periods, and the Holocene, between the early and mid-late Holocene periods. Differences in the median metric values between the periods were calculated for rarefaction, beta diversity,°disorder and skewness. To test whether the magnitude of change between the observed medians exceeded the expectations for random noise, repeat analyses were run on 10 000 randomised datasets ( Supplementary Figs. 8 and 9).
To test whether the assemblages in the later sections of the records (Younger Dryas and mid-late Holocene) were on a different trajectory to the earlier sections of the records (Bølling-Allerød and early Holocene), an ARIMA model was applied to the Bølling-Allerød and early Holocene observed metric values, using the R function arima (p, d, q), in which p was the autoregressive order, d was the degree of differencing, and q was the moving-average order (Hyndman et al., 2020). If the observed Younger Dryas and mid-late Holocene metric values fell outside the 95% confidence interval of the ARIMA forecasts, the chironomid assemblage structure had undergone substantial change between the Bølling-Allerød to Younger Dryas and the early to mid-late Holocene subsections.
To establish the relationship between the chironomid assemblage structure and climate, we assessed the correlation between the ecological metrics and the Greenland NGRIP ice core isotope record. Isotopic values were calculated for the chironomid sample dates using the predict function in R statistical software v.3.6.0 (R Core Team, 2019) ( Supplementary Fig. 10). To assess the correlation between the metric and isotope trends, Pearson's correlation coefficient was obtained from the association of the calculated isotopic values with each of the metric outcomes for the empirical datasets. To test the correlation between the metrics and isotope fluctuations, Pearson's correlation coefficient was obtained from the association between the detrended fluctuations in the ecological metrics and isotope record, using the first differences in the metric outcomes and isotopic values, obtained with the diff function in R statistical software v.3.6.0 (R Core Team, 2019). To test for non-random correlation, observed magnitudes were ranked within the set of 10 000 repeat analyses run on randomised datasets ( Supplementary Figs. 11-14).

Changes in chironomid community structure during the Late Glacial
The Late Glacial chironomid records indicated similar ecological trends, as indicated by the DCA (Fig. 3, upper panel). The large degree of variation (c. 3 σ) on axis 1 indicated large taxon turnover, from the Bølling-Allerød to the Younger Dryas. The Late Glacial chironomid assemblages plotted within the Norwegian calibration dataset (Fig. 3, lower panels), indicating that the Norwegian calibration dataset has good analogues for the fossil datasets. Goodness of fit analyses demonstrated that for the Ashik record 95.6% of samples had a good fit, for the Abernethy record 88.7% of samples had a good fit, and for Muir Park 93.0% of samples had a good fit (<95% of calibration-set squared residual distances; Supplementary Fig. 4). The Late Glacial samples plotted across the temperature contours calculated from the Norwegian calibration (Fig. 3, lower panels), suggesting that these chironomid records had been influenced by a range of temperatures over the record.
Rarefaction, beta diversity and network skewness revealed changes in assemblage structure from the Bølling-Allerød to Younger Dryas in the Late Glacial records that were largely consistent between sites (Fig. 4). Ashik and Abernethy showed decreases in rarefaction and all three sites showed increases in beta diversity and skewness from the Bølling-Allerød to Younger Dryas.°Disorder indicated a mixed response, decreasing in Ashik, increasing in Abernethy, and not changing detectably in Muir Park.
ARIMA forecasting indicated whether trends in the ecological metrics for the Younger Dryas were predictable from the Bølling-Allerød trends (Fig. 5). Rarefaction decreased during the Younger Dryas in the Ashik and Muir Park records. Beta diversity increased in all three records.°Disorder decreased during the Younger Dryas in the Ashik and Muir Park records but remained within the lower boundary of the 95% prediction interval. In the Abernethy record,°disorder increased during the Younger Dryas, remaining largely within the upper boundary of the 95% prediction interval. Skewness remained predominantly within the 95% prediction interval for all three records.
Some records indicated correlations between the ecological metrics and NGRIP isotope record trends (Fig. 6). Rarefaction at Ashik and Abernethy correlated positively with the isotope record. Beta diversity in all three records was negatively correlated with the isotope record.°Disorder did not correlate with the isotope trend in any record. Skewness in the Ashik and Abernethy records correlated negatively with the isotope record. None of the detrended datasets showed evidence of correlation between the ecological metrics and NGRIP isotope record ( Supplementary Fig. 15).
Changes in chironomid community structure during the Holocene DCA of the Holocene records are shown in Fig. 7 (upper panels). There was a large degree of variation across axis 1 (c. 2.5 σ) in the Horntjernet record, suggesting the primary driver, likely temperature, explained a large proportion of the assemblage change from the early and mid-late Holocene assemblages. The Bjornfjelltjønn and Holebudalen records showed less variation across axis 1, indicating that the primary driver was less influential on these assemblages. Some of the early Holocene samples in the Horntjernet and Bjornfjelltjønn records plotted separately from the bulk of the samples, suggesting that other drivers may have influenced the assemblages at this time. The Holocene chironomid assemblages plotted within the Norwegian calibration dataset (Fig. 7, lower panels), indicating that the Norwegian calibration dataset had good analogues for the fossil datasets. Goodness of fit analyses demonstrated that for the Horntjernet record 84.5% of samples had a good fit, for the Bjornfjelltjønn record 98.9% of samples had a good fit, and for Holebudalen 96.2% of samples had a good fit (<95% of calibration-set squared residual distances; Supplementary Fig. 5). The chironomid records plotted across a small number of the temperature contours calculated from the Norwegian calibration dataset (Fig. 7, lower panels), suggesting that the assemblages experienced a lower magnitude of temperature change than the Late Glacial assemblages. The Horntjernet record plotted across a larger number of temperature contours suggesting that temperature had a greater influence on the Horntjernet record than the Bjornfjelltjønn or Holebudalen records.
Rarefaction, beta diversity,°disorder and skewness revealed fewer detectable trends across the Holocene than during the Late Glacial (Fig. 8). Rarefaction decreased from the early to mid-late Holocene in the Horntjernet and Holebudalen records. Beta diversity decreased in the mid-late Holocene in the Horntjernet and Bjornfjelltjøn records.°Disorder increased in the Bjornfjelltjøn record. Skewness did not change in any of the records.
The ARIMA forecasting models indicate that trends in ecological metrics in the mid-late Holocene were mostly predictable from the trends in the early Holocene. The majority of the rarefaction, beta diversity, and skewness values for the mid-late Holocene fell within the prediction boundaries from the early Holocene (Fig. 9).°Disorder exceeded the predicted values in all three records. This inference from the statistical analysis is, however, influenced by the small amount of data available for the early Holocene (due to the requirement for a sample window in which to calculate°d isorder). Evidence for biometrics tracking the NGRIP isotope trend varied between metrics and locations in the Holocene (Fig. 10). Beta diversity in all three datasets negatively correlated with the NGRIP isotope trend.°Disorder correlated negatively with the NGRIP isotope trend in the Bjornfjelltjønn record. Few of the datasets indicated correlation between the detrended ecological metrics and NGRIP isotope records ( Supplementary  Fig. 16).

Discussion
Ecosystem resilience and adaptability in the face of environmental change are affected by the rate and/or magnitude of the change experienced (Overpeck et al., 1991;Skelly et al., 2007;Grimm et al., 2013). The different rates and magnitudes of the climatic changes within the Bølling-Allerød to Younger Dryas transition and Holocene were reflected in the chironomid assemblage responses.

Ecosystem structural change during rapid, high-magnitude climatic change
The chironomid assemblage compositions and structures changed in response to the abrupt, high-magnitude climate change at the Bølling-Allerød to Younger Dryas transition (Rasmussen et al., 2006;Golledge, 2010). This was evident in all four metrics: rarefaction, beta diversity,°disorder and  skewness. However, the extent of the response was variable across the metrics. There were decreases in rarefaction (Ashik and Muir Park) and increases in beta diversity (all three records) trends, unpredicted by the ARIMA models, suggesting that the changes in composition were unlikely without external forcing. Decreases in rarefaction were also detected in association with the Younger Dryas in Late Glacial chironomid records across Europe (Engels et al., 2020), further indicating that the abrupt, high-magnitude climate change experienced at the Bølling-Allerød to Younger Dryas transition was influential on chironomid communities. The correlation between the beta diversity and NGRIP ice core isotope trends suggests that taxon turnover was not independent of climate change. While it is logical that turnover should occur in response to climate change as chironomids are highly temperature sensitive (Brooks et al., 2007), these analyses confirm that high-magnitude changes in stress can drive high-magnitude responses in taxonomic composition.°D isorder and skewness were used here to indicate changes within the ecosystem structure, specifically taxonomic organisation and connectivity (Doncaster et al., 2016;Wang et al., 2019;Mayfield et al., 2020). In the Late Glacial records,°disorder displayed inconsistent trends, decreasing in Ashik, increasing in Abernethy, and indicating no clear change in Muir Park. A rise in°disorder values indicates greater unpredictability between samples, i.e. fewer taxa in common, while a decrease in°d isorder suggests greater predictability, i.e. more taxa in common (Doncaster et al., 2016). Therefore, these records indicate that changing assemblages, with decreased richness and increased turnover, can produce different outcomes in°disorder. Doncaster et al. (2016) theorised that the ordered or disordered turnover of an assemblage corresponds to the functional type of the taxa gained or lost from the assemblage. For our sites, this suggests that the functional types of the taxa lost, or gained, was different between the Ashik and Abernethy records, causing the opposing changes in°disorder values. For example, the Ashik and Abernethy records had similar assemblage compositions during the Younger Dryas, with assemblages dominated by the cold-stenothermic Micropsectra radialis-type (Chase et al., 2008;Rolland et al., 2008). However, the two lakes had different Bølling-Allerød assemblages (Brooks et al., 2012b). In the Ashik record, Microtendipes pedellus-type and Corynocera ambigua dominated the assemblage during the Bølling-Allerød, indicating relatively warm, mesotrophic conditions (Brodersen and Lindegaard, 1999;Watson et al., 2010), and the possible presence of aquatic macrophytes as suggested by the presence of Dicrotendipes nervosus-type (Brodersen et al., 2001;Watson et al., 2010). In the Abernethy record, the Bølling-Allerød assemblages were more diverse, with the continued presence of cool-adapted taxa and the arrival of more phytophilic taxa, such as Psectrocladius and Tanytarsus glabrescens-type (Brodersen et al., 2001;Heiri and Lotter, 2010). This indicates that the lakes had different lake habitats, and thus different assemblages, prior to undergoing the Younger Dryas climatic cooling, perhaps explaining the different responses within the°disorder metric. This also highlights the importance of using long-term palaeoecological records to understand the historical contingencies of communities, in order to anticipate future responses to climate and environmental change (Hawkes and Keitt, 2015).
Skewness rose in all three Late Glacial records, possibly as a result of reduced interspecific competition indicated by the reduced rarefaction (Eurich et al., 2018). A more positive skewness indicates a greater proportion of weakly connected taxa (Wang et al., 2019). It is thought that undisturbed networks are positively skewed, indicating stable, highly evolved, self-organised, complex states, where the majority of taxa are weakly connected with a small number of highly connected taxa (Albert and Barabási, 2002;Scheffer et al., 2012). This structure provides a degree of resilience to  taxonomic loss and network fragility, assuming that the weakly connected taxa are lost first (Albert et al., 2000;Sole and Montoya, 2001;Dunne et al., 2002). Wang et al. (2019) indicated a reduction in positive skewness in association with increased environmental stress, thus here we hypothesised that an increase in climatic stress would simulate a decrease in skewness. In the records studied here, organic matter content and head count concentration decreased during the Younger Dryas (Brooks et al., 2012b(Brooks et al., , 2016, suggesting lower biological productivity and more homogeneous environments. The decrease in taxon richness (rarefaction) may have reduced the competition for the remaining micro-niches, thus reducing the number of taxon interactions and causing the rise in skewness.
The rise in positive skewness during the Younger Dryas was relatively small, remaining within the prediction intervals of the ARIMA model. This suggests that, despite the large compositional turnover, the change in ecosystem structure was relatively small. This may relate to the ecological function of the lost taxa within the ecosystem. Gallagher et al. (2013) found that loss of taxa does not necessarily lead to impairment of ecosystem function as the remaining taxa may still occupy all the available niches, while Obertegger and Flaim (2018) found that functionally similar ecosystems can exist in large and small communities, at different levels of diversity, based on resource use. Furthermore, many chironomid taxa have relatively large trait variability or plasticity (Serra et al., 2017), and depending on availability of food resources, taxa can display generalist or opportunistic feeding behaviours (Reuss et al., 2013;Lee et al., 2018), thus increasing the ability of different chironomid taxa to substitute functions within an ecosystem and reducing the effect of taxon loss on the ecosystem functionality (Brown et al., 2001). At genus level, there is an overlap in food and habitat preferences within Orthocladiinae and Chironominae (sub-tribes Chironomini, Tanytarsini and Pseudochironomini) (Serra et al., 2016). The assemblages at Loch Ashik, Abernethy and Muir Park primarily consisted of these Chironomidae genus-types, thus the temperature-driven assemblage changes may have had a subdued effect on ecosystem functionality.
The majority of the chironomid samples in this study produced negative skewness values, suggesting that the communities were prevented from developing a more complex structure. The lack of positive skewness in chironomid communities could relate to unstable environmental conditions, or non-analogue communities when comparing past conditions with modern assemblages (Velle et al., 2005b). However, the lack of positive skewness could also be a limitation of the analytical method as the absolute skewness values are sensitive to statistical procedures. Mayfield et al. (2020) also reported a predominance of negatively skewed chironomid communities in their study of structural change in chironomid communities across high-latitude regions; thus here we focused on the relative trends as a guide to structural change.
Additional drivers influence ecosystem structure during periods of low climatic change In the Holocene, the change in beta diversity was relatively gradual, reflecting the lower-magnitude changes in the NGRIP record (Wanner et al., 2008). There was a rise in°disorder in the Bjornfjelltjønn record indicating a rise in composition dissimilarity with the changing climate. Skewness indicated little ecosystem structural stability. Theoretically, skewness should become more positive during long, stable periods, such as the Holocene, as habitats develop and micro-niches diversify (Wang et al., 2019). The lack of clear trends in the Holocene records is probably an effect of non-climate factors, as discussed below.
As the climate stabilised in the mid-late Holocene (Wolff et al., 2010), other environmental drivers may have become prominent, such as edaphic factors, paludification, vegetation change, and more recently, anthropogenic drivers including increased nutrient loading (Velle et al., 2010). Such drivers can also cause changes in ecosystem structure and functionality; for example, through eutrophication and acidification (Heiri and Lotter, 2003;Ilyashuk et al., 2005;Velle et al., 2010). While ecosystem stability may be ultimately controlled by climate (Birks and Ammann, 2000;Brooks et al., 2012b), our study suggests that localised factors may have a greater influence on assemblage structure in times of low climatic stress, as also suggested by Velle et al. (2005a) and Brooks (2006). Inspection of the chironomid assemblages supports this: Horntjernet, situated at the current boreal ecotone boundary in arctic Norway, appears to have been exposed to post-glacial environments with an early dominance of coldstenothermic taxa, followed by taxa suggesting subsequent paludification and afforestation during the mid-Holocene (unpublished data). Landscape change and increased dissolved organic carbon possibly caused the large change in°d isorder at c. 5000 to 2000 cal a BP. Bjornfjelltjønn, a lowalpine lake in northwest Norway, was subjected to increased meltwater during the early Holocene (suggested by the presence of Abiskomyia) and acidification during the later Holocene in response to vegetational succession and soil development (indicated by the presence of Heterotrissocladius, Zalutschia and Heterotanytarsus taxa types) (Brooks, 2006). Plant macrofossils indicated a change in landscape vegetation at c. 3800 cal a BP (Brooks, 2006), possibly causing the brief decrease in°disorder in the chironomid record around this time. Holebudalen, a lowalpine lake from southern Norway, also experienced acidification during the Holocene, most likely due to local forest or bog development (indicated by the presence of Zalutschia, Psectrocladius sordidellus-type and Sergentia coracina-type) (Velle et al., 2005a;Brooks et al., 2012a). Velle et al. (2005a) identify a shift in chironomid composition at c. 2000 cal a BP, attributing this change to an increase in thermophilic and eutrophic taxa. Thus, a change in the chironomid assemblage to taxa more sensitive to water quality may have caused the changes in the°disorder metric.
Other studies focusing on Holocene chironomid records have indicated complex relationships between chironomid community diversity and changes in climate and/or environment, e.g. Velle et al. (2005a) and Engels et al. (2020). Comparisons with proxies of landscape dynamics and water quality, such as pollen or diatoms, could help clarify drivers of assemblage change during periods of lower climate change and provide greater insight into changing lake ecosystem structure; for example, as done by Rosén et al. (2001), Langdon et al. (2004) and Velle et al. (2005b). This emphasises the importance of constructing detailed, multiproxy analyses of lake records to enable full comprehension of the changes preserved in the records. Furthermore, the low-amplitude climate change experienced during the Holocene was spatiotemporally heterogeneous across Europe (Wanner et al., 2008;Renssen et al., 2009), and thus the lakes studied here, located on a latitudinal gradient, may have experienced different rates and magnitudes of climate change. Comparison with Holocene-specific climate reconstructions, such as Kaufman et al., (2020a) and (2020b), may provide greater insight into the effect of localised Holocene climate dynamics.

Timescales and magnitude of change
The selected proxy and sampling resolution must have the appropriate sensitivity to record the magnitude of change and the response of a system to a driver of stress (Lenton, 2011;Lenton et al., 2012). The records used in this study span millennial timescales, with variable temporal resolution; the time between adjacent samples represented c. 3-163 years during the Bølling-Allerød and c. 43-369 years in the Younger Dryas.
Similarly, in the Holocene, the temporal resolution varied between c. 44 and 415 years in the early Holocene and c. 36-430 years in the mid-late Holocene. While these temporal resolutions were comparable, changes in°disorder and skewness were identified in the Late Glacial records, but not in the Holocene records. This indicates that rapid, high-magnitude climate change can drive ecosystem structural change, as seen at the Bølling-Allerød to Younger Dryas transition. However, the slower, lower-magnitude, climate change between the early and mid-late Holocene was not great enough to influence the ecosystem structure. Furthermore, the metrics did not indicate substantial changes in the trends in association with rapid, higher-magnitude Holocene climatic events, such as the Preboreal Oscillation, 9.3 ka and 8.2 ka events (Alley et al., 1997;Rasmussen et al., 2007) Fig. 17). Previous studies have indicated changes in chironomid assemblage composition and associated chironomid-inferred temperature reconstructions with these Holocene climatic events (Axford et al., 2008;Larocque-Tobler et al., 2010;Paus et al., 2011;Porinchu et al., 2019). This suggests that while chironomid assemblages can be sensitive enough to record taxonomic changes in association with such events, it is possible that these resulting assemblage changes had a limited effect on the ecosystem structure, or community structural changes were not detectable in the temporal resolutions presented here. Other environmental stressors have been seen to drive structural change on shorter timescales; Doncaster et al. (2016) and Wang et al. (2019) found high-magnitude changes in°disorder (chironomids and diatoms) and skewness (diatoms) in lakes experiencing high-magnitude shifts in total phosphorus that drove large changes in chironomid and diatom assemblage structure on a sub-decadal timescale. This reaffirms that the sampling resolution must have the appropriate sensitivity to record the magnitude of change and the response of a system to a driver of stress (Lenton, 2011;Lenton et al., 2012). Analysing records with a higher temporal resolution spanning rapid climatic events could provide greater insight into the ability of rapid climate events to drive structural change.

Future possibilities of ecosystem structural analysis
During the Late Glacial, the climate experienced a rapid, high-magnitude decrease in temperature (Dansgaard et al., 1993;Wolff et al., 2010), thus this paper quantifies the effects of abrupt, high-magnitude climatic cooling on chironomid community structure. Considering that the current climate is warming (Smith et al., 2015), a further study testing chironomid structural change from cool to warm climates may be informative; such as the Younger Dryas to Holocene transition. This was not done as part of this study due to a lack of sample availability; the Younger Dryas subsections comprised a small number of samples (Ashik and Abernethy n = 7, Muir Park n = 8). It was thought that these records do not contain enough samples to generate a reliable comparison between the Younger Dryas and Holocene, particularly considering that the calculation of°d isorder requires a window of samples. Processing and analysing cores with greater sampling resolution during the Younger Dryas to Holocene transition may provide the opportunity to quantify structural change from cool to warm climatic conditions and increase our understanding of the effect of rapid climatic warming on chironomid community structure.
The records analysed here provided insight into the influence of past millennial-scale climatic changes on the structure of chironomid communities, yet to fully understand how ecosystem structures may respond to future rates and magnitudes of climate change, we need to consider how ecosystems have responded to stress in the more recent past. Chironomid communities have already experienced changes in composition in relation to recent 20th century warming (Porinchu et al., 2007); however, finding lakes without additional environmental or human-derived stress is increas-ingly unlikely, even in the Arctic (Smol et al., 2005). Considering that the Holocene records studied here are likely to be influenced by multiple interacting factors, this further emphasises the importance of investigating the effect of combined environmental stressors. We propose that analysing high-resolution lake records spanning the more recent past may further our understanding of the combination of ecosystem drivers that lakes are currently facing and the effects such drivers have on lake ecosystem structures; for example, Ilyashuk et al. (2015), Nevalainen et al. (2015) and Engels et al. (2020).
Macroinvertebrates, such as chironomids, are a fundamental component of freshwater ecosystems (Jones and Grey, 2004;Ólafsson and Paterson, 2004), and changes in the abundance or distribution of the invertebrate component of a food web can have repercussions throughout the ecosystem (Petchey et al., 1999). Multi-trophic scale ecosystems are at risk from climate change (Beier, 2004;Pecl et al., 2017), and analysis of the interactions between trophic levels has indicated that low trophic levels can indicate the level of deterioration in food-web stability (Kuiper et al., 2015;Kivilä et al., 2019). It is thought that Figure 10. Correlation between the NGRIP isotope trends (black) and the chironomid ecological metrics (red) for the Holocene records. Calculated isotopic values are shown for the same ages as the chironomid samples. Pearson's correlation coefficients indicate the correlated trend between each metric and the isotopic values. Asterisks denote p < 0.05, where the observed magnitude of correlation was in the top 5% of correlations obtained from 10 000 replicate datasets with a randomised order of quantities. [Color figure can be viewed at http://wileyonlinelibrary.com] chironomids can improve biomonitoring approaches and increase understanding of lake ecosystem change (Nicacio et al., 2015;Czechowski et al., 2020). It is hoped that through improving our knowledge of chironomid community sensitivity to environmental change with this study, we can contribute to the wider understanding of lake ecosystem resilience and stability.

Conclusions
This study evaluated the effect of the rate and magnitude of climate change on chironomid community structure using two periods of past climatic change; the relatively rapid, highmagnitude climate change at the Bølling-Allerød to Younger Dryas transition, and the more gradual, relatively lowmagnitude climate change during the early-mid-Holocene transition. The Late Glacial records exhibited high-magnitude changes in beta diversity, but changes in the community structure, represented by°disorder and skewness, were smaller. This was perhaps due to the functional resilience of the chironomid community through the replacement of same functional-type taxa within the network. The Holocene records indicated low taxonomic turnover in association with low levels of climate change, with little indication of climatedriven structural change during the Holocene. It is likely that other environmental factors were having a greater influence on the chironomid assemblage at this time. Overall, this study illustrates that beta diversity is sensitive to compositional change driven by high-and low-magnitude temperature change. High-magnitude changes are required to cause structural change within chironomid communities, as seen in the Late Glacial records. The lack of changes in°disorder and skewness in the Holocene records suggest that low-magnitude climate change is not a key driver of structural change. We emphasise that a greater understanding of how combined environmental stresses may influence community structure is important for understanding current and future lake ecosystem change.
Author contributions-R. J. Mayfield, P. G. Langdon, C. P. Doncaster and J. A. Dearing discussed the research conceptualisation and outcomes. C. P. Doncaster provided the original R code for the compositional disorder calculations. R. Wang provided the original MATLAB network skewness code. R. J. Mayfield adapted the above codes and ran all analyses on the empirical data. R. J. Mayfield was part of the coring team at Horntjernet, and processed and identified the Horntjernet chironomid samples at the University of Southampton. P. G. Langdon and S. J. Brooks provided assistance with chironomid identification. Other chironomid data were provided by S. J. Brooks, K. Davies, and G. Velle. R. J. Mayfield wrote the draft manuscript, on which all co-authors commented.

Supporting information
Additional supporting information can be found in the online version of this article.
Acknowledgements. This study was supported by a PhD studentship awarded to R. J. Mayfield provided by the UK National Environmental Research Council (grant no. NE/L002531/1). We would like to thank I. G. Alsos for organizing fieldwork at Horntjernet, D. Rijal and F. J. A. Murguzur for their fieldwork assistance, and P. Heintzman and D. Rijal for producing the Horntjernet age model. The fieldwork was supported by the project 'ECOGEN -Ecosystem change and species persistence over time: a genome-based approach' at The Arctic University Museum of Norway (Research Council of Norway grant number 250963/F20 to I.G. Alsos). We would like to extend our thanks to all contributors of the chironomid datasets, including C. T. Langdon and G. Schellinger who provided assistance with chironomid identifications for the Horntjernet record. Thank you to A. Dugmore for discussions on calculating rates of change. We thank Stefan Engels and a second anonymous reviewer for their valuable comments and helping us improve this paper.

Data availability statement
Chironomid data are available from the original authors, see Table 1. The beta diversity code is available from the R documentation for the adespatial package (Dray et al., 2019). The network skewness MATLAB code is available online in the supplementary information from the original publication (Wang et al., 2019). R code for calculating compositional disorder can be found in the supplementary information for this paper.