Stability of chironomid community structure during historic climatic and environmental change in subarctic Alaska

By understanding lake ecosystem resilience in the face of increasing environmental and anthropogenic stress, we can hope to anticipate future ecosystem instability. We assess recent historic ecosystem resilience using composition and network analyses of empirical zoobenthos chironomid (Diptera: Chironomidae; nonbiting midges) reconstructions from three Subarctic Alaskan lakes, spanning the last c. 200 yr. We measured community richness, turnover and structure using taxon richness, beta diversity, and network skewness, respectively. Simulated taxonomic networks were created to establish the sensitivity of these metrics to changes in taxon connectivity, and to inform the interpretation of empirical chironomid records. The models indicated that beta diversity was more sensitive to taxon loss, while skewness was more sensitive to taxon gain. Both beta diversity and skewness were required to understand structural change under taxon replacement. The simulated arrival of strongly connected taxa caused a greater decrease in skewness than the arrival of weakly connected taxa. The empirical data sets indicated a rise in taxon richness (measured as rarefaction) and beta diversity in the recent samples. Changes in chironomid composition were associated with climate warming (replacement of cold taxa with temperate taxa) and increased lake biological productivity (the arrival of macrophyte‐associated taxa). Skewness was predominantly negative across the lakes, indicating high taxon connectivity and structural stress. However, little directional change in the skewness trends suggests some resilience within the chironomid community structures in relation to the current levels of climate and environmental stress. Continued climatic warming, and associated rises in nutrient levels, may cause further structural stress and ecological degradation.

The world is warming at an unprecedented rate (Smith et al. 2015), with high-latitude regions warming almost twice as fast as the global average (Walsh et al. 2017;IPCC 2018). High-latitude freshwater ecosystems are particularly vulnerable to climate change, with recent climate-driven changes in lake assemblages already recorded in Arctic and Subarctic locations (Smol and Douglas 2007a;Rosén et al. 2009). Combined with growing anthropogenic stresses Vilmi et al. 2017), such as increased nutrient influx (e.g., from sewage inputs), high-latitude lakes have also experienced a rise in nutrient levels (Antoniades et al. 2011;Gallant et al. 2020). High-latitude lakes are important resources for local communities, providing drinking water and fisheries (Sjolander 2011), thus changes in lake ecosystem state have both ecological (Post et al. 2009;Woelders et al. 2018) and human (Hayden et al. 2017) consequences. Climate and anthropogenic stresses are known to impact ecosystem richness and turnover, but less is known about how such stressors can impact ecosystem structures. Here, we aim to improve our understanding of current freshwater aquatic ecosystem stability, under increasing climatic and anthropogenic stress, through compositional and network analyses of microfaunal data sets in lake sediments spanning the last c. 200 yr.
Climate-driven increases in the abundance and diversity of chironomid taxa have been recorded across the Arctic over the last c. 200 yr (Quinlan et al. 2005;Luoto et al. 2014). Warmer temperatures can increase lake microhabitat heterogeneity and biological productivity; for example, through increased macrophyte presence, which in turn can increase microhabitat and faunal diversity (Langdon et al. 2010;Lau et al. 2020). Initial rises in chironomid taxon richness with increased temperatures have been observed across subarctic regions (Mayfield et al. 2020). However, combined effects of climate and other environmental pressures, such as bird habitats, can have detrimental effects on chironomid distribution (Luoto and Ojala 2018;Mayfield et al. 2020). Changes in the abundance or distribution of aquatic invertebrates can have repercussions throughout the ecosystem (Quinlan et al. 2005; Pearce-Higgins 2010), such as altering ecosystem function and corresponding food web structures (Chapin et al. 2000;Wrona et al. 2006). Climate can indirectly influence assemblage structure through changing taxonomic interspecific interactions (Harley 2011), which can alter the resilience of an ecosystem Grimm et al. 2013). Here, we define ecosystem resilience as the capacity of an ecosystem to absorb disturbance and retain the same ecological function, without shifting into an alternative state (Holling 1973;Folke et al. 2004;Walker et al. 2004). In our study, we use temporal beta diversity to quantify taxonomic turnover (Legendre and Gauthier 2014) and ecosystem organization (Soininen et al. 2007;Heino et al. 2012). Community turnover under exogenous stress is often influenced by the loss of sensitive species (Pound et al. 2018). Considering the large number of sensitive, coldadapted chironomid taxa within the Arctic (Medeiros et al. 2011;Eggermont and Heiri 2012), beta diversity should effectively indicate climate-driven changes in community composition.
Biotic communities comprise networks of inter-connected taxa (Bruder et al. 2019). Under low stress conditions, ecosystems are often spatially heterogeneous, with the majority of taxa self-organizing into their preferred microhabitats ( Van de Meutter et al. 2005;Kouamé et al. 2011), reducing the number of interactions or co-occurrences between taxa, that is, these taxa are weakly connected (Dunne et al. 2002;van Nes and Scheffer 2005). A few taxa may have a large number of interactions or co-occurrences, that is, they are strongly connected (Strogatz 2001;Doncaster et al. 2016;Wang et al. 2019). This produces a network of predominantly weakly connected taxa. Networks with a high proportion of weakly connected taxa are more stable as they have a greater capacity to change gradually to external stressors, as individual taxa will change (leave or join the network) as conditions vary (Dunne et al. 2002;van Nes and Scheffer 2005;Scheffer et al. 2012). Substantial levels of stress can cause a loss of microhabitats and environmental homogenization, increasing the number of interactions between taxa. This reorganization of taxa causes an increase in the proportion of strongly connected taxa. In networks of strongly connected taxa, the fate of each individual taxon depends on the fates of its many other interconnected taxa. Although high connectivity can provide a local and temporary resilience to stress, it can also indicate instability because systems dominated by strongly connected taxa have an increased likelihood of systemic critical transition: a sudden and systemwide shift in taxa (Dunne et al. 2002, van Nes and Scheffer 2005, Scheffer et al. 2012). Thus, measuring community structure, that is, the distribution of taxa and their connections, can indicate the level of stress a community might be experiencing and the likelihood for transitions (Wang et al. 2019).
In this paper, we quantify the current ecosystem stability of temperature-sensitive and environmentally sensitive chironomid communities from three lakes in Subarctic Alaska, an area that has experienced climate warming (Fig. 1A,B) (Hinzman et al. 2005) and human population growth ( Fig. 1C) (World Population Review 2020) over the last c. 100 yr. Taxon richness, beta diversity, and network skewness were used to characterize ecosystem richness, taxonomic turnover, and assemblage structure in empirical chironomid records spanning the last c. 200 yr. To determine the potential for ecological July, and August) across Alaska (gray lines, see Supporting Information Fig. S1 for temperature site locations), with mean (red) and smoothed regression (black) curves. Temperature data are from the CRU TS4 Google Earth interface (Harris et al. 2020). (B) To indicate the rate and magnitude of climate change, the cumulative deviation from the long-term mean was calculated for each temperature location (gray lines), following Dugmore et al. (2007). The long-term mean was calculated over periods of 10 yr (10 samples). Mean cumulative deviations (red) and smoothed regression (black) curves are also shown. (C) Human population data are shown for the whole of Alaska (black dashed line) and the Fairbanks North Star Borough (solid red line). Population data are from the World Population Review (2020). metrics to inform on changes in taxon connectivity, and thus community structure, the metrics were also applied to simulated taxonomic networks with changing taxonomic connectivity. Trends identified in the simulated data sets will inform assessment of changes in the empirical chironomid data sets under recent climate change, indicating whether the empirical chironomid communities were showing signs of structural stress. For each lake, sedimentological and geochemical records were used to identify changes in the lake catchment and in-lake processes. Temperature observations covering the last c. 100 yr were used as an indicator of recent climate change. Ordination was used to assess the influence of environmental variables on the chironomid assemblages. We hypothesize that if the chironomid communities have been affected by recent climatic warming, taxonomic richness and beta diversity (turnover) are likely to have increased. Network skewness should indicate whether this warming had a positive or negative effect on the community structure. For example, if the chironomid communities were experiencing temperature stress, we may expect to see a loss of weakly connected taxa expressed by a shift from positive toward negative skewness values as the most vulnerable (cold stenothermic) taxa disappear from the assemblages.

Changing taxon connectivity in hypothetical networks
Increased environmental stress can increase taxonomic connectivity (i.e., the number of connections a taxon has) as vulnerable species become extinct and/or new taxa colonize, resulting in structural reorganization (Wang et al. 2019). To assess the informative potential of taxon richness, beta diversity and network skewness in detecting changes in taxonomic connectivity, the metrics are applied to simulated taxonomic networks with changing taxonomic connectivity. Connectivity is defined as the associations of taxa co-occurring in the same microhabitat (Wang et al. 2019;Mayfield et al. 2020). Hypothetical models have previously been used to produce null expectations of community composition and structural change across temperature gradients in Mayfield et al. (2020); however, those models did not account for the connectivity between taxa. Here, we test the effect of changing connectivity on the metrics.
Networks with a scale-free distribution, that is, a large proportion of weakly connected taxa and a small proportion of strongly connected taxa (Barab asi and Bonabeau 2003), were generated to represent a single lake community under low stress conditions. Networks were generated using the barabasi.game function from the igraph package (Csardi and Nepusz 2006) for R statistical software v. 3.6.0 (R Core Team 2019). Taxa correspond to a species morphotype and connections refer to undirected, unweighted taxon co-occurrences within lake microhabitats. Changing taxon presence is likely to occur as environmental conditions change under exogenous stress and taxa leave or join the community network. However, the type of environmental stress driving the taxon change is not discussed in these model scenarios. The effects of changing taxon connectivity under taxon loss, gain, and replacement were tested using scenarios of decreasing and increasing taxonomic richness and taxon turnover (Table 1). At each iteration, the model randomly selected the number of taxa (1 or 2) to change within the network. Starting network graphs of 100 taxa were generated for the taxon loss (scenarios 1-3) and replacement (scenarios 7-11) scenarios. Starting graphs of 10 taxa were generated for the taxon gain scenarios (scenarios 4-6). Each model had 50 iterations, with each iteration representing a time step. Fifty iterations were selected to show the general trends of changing taxon presence and maintain a manageable network size. The output was a presence-absence matrix indicating the changing presence of taxa, and their number of connections (called taxon degree), within the network over time. Taxa with zero connections were removed from the network, ensuring that all taxa had at least one connection to another. The network could become fragmented or additional taxa could be lost as a consequence of the removal of a taxon with multiple connections. Each scenario was repeated 100 times to create 100 output matrices. All analyses were run in R Studio v. 3.6.0 (R Core Team 2019).
In scale-free networks, the distribution of taxon degrees follows a power law (Barab asi and Albert 1999). Under linear preferential attachment, taxa joining a network have a greater likelihood of connecting with strongly connected taxa. A power law of 1, for linear preferential attachment, was used to calculate the new degree distribution as taxa were added in scenarios 4, 5, and 7-10. In scenarios 6 and 11, a power law of 0 was used to enable newly arrived taxa to join a random node within the network. Taxa were removed from the existing networks using the delete_vertices function and added to the networks using the sample_pa function from the igraph package. The sample_pa function requires an algorithm to generate the new network structure. The "bag" algorithm was used for scenarios 4, 5, and 7-10; "bag" works by placing the applicable taxa into a multiset the same number of times as their (in-) degree plus one, and by drawing the required number of cited vertices from the bag, with replacement (Csardi and Nepusz 2006). Under this method, multiple edges between the same two taxa could be generated. Multiple connections between taxa were removed from the networks using simplify. The "bag" algorithm cannot work with a power law of 0, thus, the "psumtree" algorithm was used for scenarios 6 and 11. The "psumtree" algorithm uses a partial prefix-sum tree to generate the new network structure (Csardi and Nepusz 2006).
Taxon richness, beta diversity and network skewness were calculated on the simulated output matrices. Taxonomic richness was calculated as the sum of taxa present in each network iteration. Beta diversity was calculated using the variance partitioning framework by Legendre and De Caceres (2013), using the beta.div function and "Hellinger" method in the adespatial package in R (Dray et al. 2019). To measure skewness, first the frequency distribution of connections between the taxa in the network was calculated, and then the skewness of the frequency distribution was calculated using the skew function, from the psych package for R (Revelle 2016).

Chironomid community change in three Alaskan lakes Site selection and descriptions
Three lakes from Subarctic Alaska were sampled and analyzed to ascertain the current ecosystem stability of environmentally sensitive chironomid communities ( Fig. 2A; Supporting Information Table S1). The lakes are part of the North American chironomid calibration data set (Fortin et al. 2015) and have previously been sampled by Barley et al. (2006). The lakes were originally selected for their small size, shallowness, and lack of disturbance, including little or no inflow, to ensure the chironomid assemblage reflected local environmental and limnological conditions (Barley et al. 2006). Lakes for this study were selected based on location, accessibility, lack of disturbance, and comparable environmental conditions. The lakes are located within the boreal and discontinuous permafrost zones of Alaska (Brandt 2009;Wang et al. 2018) and situated on unconsolidated Quaternary surficial deposits . The lakes have experienced climatic warming ( Fig. 2B) with comparable rates of change (Supporting Information Fig. S2) and trends during the last c. 100 yr. The differences between the microclimates may relate to altitude; Round Tangle is located at the highest altitude (877 a. s.l.), while Little Lost lake has the lowest altitude (303 a.s.l.) (Supporting Information Table S1).
Round Tangle lake is part of a chain of lakes (c. 26 km long) connected by streams. Samples were taken from the smaller (c. 0.04 km 2 ), shallower southwest embayment, located by the campground. The surrounding catchment was largely comprised of shrub tundra, dominated by dwarf birch (Betula glandulosa), various willow (Salix) species, and scattered white spruce (Picea glauca).
East Cobb lake is a shallow, c. 1.1 km 2 , lake situated just off the Glen Highway, near the small settlement of Slana (population: c. 150 as of the 2010 census, World Population Review 2020). The catchment vegetation is dominated by white spruce (Pi. glauca), Alaska birch (Betula neoalaskana), green alder (Alnus viridis), balsam poplar (Populus balsamifera), and various willow (Salix) species. The lake edges were lined with grasses. Pondweeds (Potamogeton) and mollusk shells were visible in the shallower sections of the lake.
Little Lost lake is a small, c. 0.3 km 2 , fishing lake, next to the larger Quartz Lake, situated in the Quartz Lake State Recreation Area, c. 10 miles north of Delta Junction. The catchment vegetation is dominated by white spruce (Pi. glauca), with some Alaska birch (B. neoalaskana), gray alder (Alnus incana), and various willow (Salix) species. Yellow (Nuphar lutea) and white (Nymphaea alba) water lilies were present within the lake.

Core collection and analysis
Surface sediment cores were collected from the lakes, summer 2018, using a percussion gravity corer and polycarbonate tubes. A single-beam echo sounder was used to ascertain the Table 1. Scenarios for testing the effect of taxon connectivity and taxonomic loss, gain, and replacement on three ecological metrics: taxon richness, beta diversity, and network skewness. For the taxon replacement scenarios, taxa were removed or added to the network at alternate iterations.

Taxon loss
Scenario 1: Weakly connected (1-2 connections) taxa were removed from the network Scenario 2: Strongly connected taxa were removed from the network. The three most strongly connected taxa could be selected for removal Scenario 3: Taxa were removed from the network, irrespective of the number of connections Taxon gain Scenario 4: Weakly connected (1-2 connections) taxa were added to the network Scenario 5: Strongly connected taxa were added to the network. A newly added taxon could have the same number of connections as one of the 3 most strongly connected taxa within the existing network. This meant that the number of connections the new taxa had was proportional to the existing network Scenario 6: Taxa were added to the network with a random number of connections (1-24 connections) Taxon replacement Scenario 7: Weakly connected (1-2 connections) taxa were replaced with other weakly connected (1-2 connections) taxa Scenario 8: Strongly connected taxa were replaced with other strongly connected taxa. The three mostconnected taxa could be selected for removal, while taxa added to the network could have the same number of connections as the remaining three most strongly connected taxa Scenario 9: Weakly connected (1-2 connections) taxa were replaced with strongly connected taxa. The taxa added to the network could have the same number of connections as the 3 most strongly connected taxa Scenario 10: Strongly connected taxa were replaced with weakly connected (1-2 connections) taxa. Removed taxa were randomly selected from the three most strongly connected taxa Scenario 11: Taxa were randomly removed and added to the network, irrespective of number of connections deepest part of the lakes. Two parallel cores were collected from each lake at the deepest part; Round Tangle 3.7 m, East Cobb 1.6 m, and Little Lost 2.7 m (see Supporting Information Table S1). For each site, one of the cores was processed for 210 Pb dating (alpha spectrometry was used to aid dating precision; however, this process is destructive causing a loss of sediment), loss on ignition (LOI) and bulk density. The second core was processed for chironomid assemblage composition, LOI, bulk density, magnetic susceptibility, and XRF sediment analysis. LOI and bulk density were used to correlate the parallel cores and align the 210 Pb dating profile with the chironomid assemblage composition (see Supporting Information Figs. S3-S5). Highresolution monthly temperature observations are available for the period 1901-2018 (Harris et al. 2020) and indicate climatic warming ( Fig. 1); however, sediment and chironomid analyses were prepared to cover the last c. 200 yr to provide increased knowledge of the abiotic and biotic lake conditions prior to the observed climate change. 210 Pb dating was used to establish a chronology for the lake cores. Analysis was carried out at GAU-Radioanalytical, Ocean and Earth Science, University of Southampton. Freeze-dried sediment samples were treated using acid digestion (Aqua Regia) with added 209 Po as a recovery tracer. The extracted 210 Po/ 209 Po was plated onto silver disks using auto electrodeposition from a diluted HCl solution. The alpha spectrometric sources were counted using ORTEC Octete spectrometers equipped with PIPS semiconductor detectors. Spectra were analyzed using Maestro software. Raw 210 Pb dates were calculated using the Constant Rate of Supply Model (Appleby and Oldfield 1978). This assumes a constant rate of supply of unsupported 210 Pb to the lake sediments. See Supporting Information Figs. S6-S8 for further information and the 210 Pb dating output.
To understand the influence of catchment dynamics and in-lake processes on the chironomid assemblage, the minerogenic and geochemical compositions of the lake sediment cores were analyzed using magnetic susceptibility and elemental X-ray fluorescence (XRF) spectroscopy at the British Ocean Sediment Core Facility (BOSCORF), National Oceanography Centre, Southampton. Magnetic susceptibility measurements were taken at 0.5 cm resolution using a Multi-Sensor Core Logger-Standard. Geochemical analysis was measured at 200 μm resolution using an ITRAX energy-dispersive XRF (ED-XRF) core scanner (Croudace et al. 2006). To account for changes in core density, surface topography, and water content, the elemental data (counts per second, cps) were normalized using the total kilo counts per second (kcps) (Croudace and Rothwell 2015): normalized elemental data = cps/kcps. LOI and bulk density analyses were carried out in the School of Geography and Environmental Science laboratories, University of Southampton. LOI was used to indicate the sediment organic matter content ) and bulk density was used to measure the compactness of the sediment matrix ( mer ambient temperatures (June, July, and August) were plotted for the 0.5 Â 0.5 latitude-longitude grids in which the three lakes are located (B). Temperature data were downloaded from the CRU TS4 Google Earth interface (Harris et al. 2020). Lake images (C-E).
density (g/cm 3 ), 0.5 cm 3 (v) samples were subsampled at a 0.5 cm resolution and dried over night at 105 C and weighed for the dry sediment mass (α). The samples were incinerated at 550 C for 2 h and weighed for the ash sediment mass (β) (Howard and Howard, 1990). Organic carbon content was calculated using the equation: LOI % = ([α -β]/α) * 100. Bulk density was calculated as the dry sediment mass (α) divided by the sample volume (v): bulk density = α/v. For magnetic susceptibility, LOI, bulk density, and ITRAX records, see Supporting Information Figs. S9-S11. Chironomid assemblage reconstructions were carried out at the School of Geography and Environmental Science laboratories, University of Southampton. Sediment samples were taken at 0.5 cm resolution and deflocculated using 10% KOH for 18 min at c. 70 C and passed through 180-and 90-μm sieves. Chironomid head capsules were picked individually from the residues using Bogorov sorting trays and a mounted pin, airdried, and mounted in Euparal. The specimens were identified to genus or species-morphotype level using standardized subfossil taxonomy (Brooks et al. 2007). For Round Tangle and East Cobb Lake, all samples had > 75 head capsules, with the exception of 1 sample (6.0-6.5 cm) from East Cobb with 71.5 head capsules (Supporting Information Table S2). For Little Lost, all lake samples had a minimum of 50 head capsules (lowest count 53.5) Quinlan and Smol 2001). Chironomid assemblage diagrams were created in R using the rioja package (Juggins 2017). Shifts in the chironomid assemblages were identified from dissimilarity indices, using the vegdist function in Vegan (Oksanen et al. 2013). Dendrographs were added to the assemblage diagrams to indicate zones of differing taxonomic composition.

Calculating community change
To identify key gradients within the chironomid assemblages, the chironomid records were plotted in ecological space defined by detrended correspondence analysis (DCA). Key environmental influences on the chironomid assemblages were identified with canonical correspondence analysis (CCA). Ordinations were run using the vegan package (Oksanen et al. 2013). Rarefaction was used as a measure of taxon richness to account for variations in count sums when comparing samples. Rarefaction was calculated using the rarefy function in vegan (Oksanen et al. 2013). Taxa counts were rounded to integers and the minimum count sum was used to calculate rarefaction for each site. Beta diversity was calculated using the beta.div function and "Hellinger" method in the adespatial package (Dray et al. 2019). The assemblage data were transformed by the "Hellinger" method within the beta.div function. The above analyses were run in R Studio v. 3.6.0 (R Core Team 2019).
Network skewness was calculated in MATLAB (ver. R2017b) following Wang et al. (2019). To calculate skewness values for temporal assemblages, a modern calibration data set is required to identify the most-frequently co-occurring chironomid taxa pairs occupying the upper two quartiles of positive values for Cramér's association coefficient (Q2 of V+). Then, taxon degree (number of connections) was calculated for the fossilized data sets based on the co-occurrences of taxa in the calibration set. Network skewness was measured as the skewness of the frequency distribution of taxon degree. Here, skewness was calculated for the fossilized data sets using both the North American (Fortin et al. 2015) and Norwegian Birks 2001, 2004) chironomid calibration data sets, as the data sets have different degrees of taxonomic resolution. To test the appropriateness of the calibration sets, nonmetric multidimensional scaling was used to ordinate the temporal chironomid records and calibration data sets simultaneously, using the Timetrack function in R (Simpson and Oksanen 2019) (Supporting Information Fig. S12). The fossil samples plotted within the calibration data sets, suggesting that both calibration data sets had good analogues for the fossil records. Goodness of fit analyses, performed using squared residual distance (Juggins and Birks 2012), also indicated that both calibration data sets were appropriate (Supporting Information Fig. S13). As required by this method of skewness, the temporal chironomid assemblages were filtered to match the taxa identified in the calibration sets, with some subcategories of taxa being grouped to genus-type. Rarefaction and beta diversity were calculated on the filtered data sets for comparison. Pearson's correlation indicated strong correlation between the rarefaction and the beta diversity trends for the full and filtered data sets. Pearson's correlation did not indicate detectable correlations between the skewness trends for the filtered data sets (Supporting Information Fig. S14).
Differences in the median metric values were calculated for rarefaction, beta diversity, and skewness for the assemblage zones (identified by the assemblage dendrographs). To test whether the magnitude of change between the observed medians exceeded the expectations for random noise, repeat analyses were run on 10,000 randomized data sets, following Telford and Birks (2011) (Supporting Information Fig. S15).

Environmental drivers of ecosystem change
Relationships between key environmental variables and the chironomid ecological metrics (rarefaction, beta diversity, and skewness) were quantified with generalized additive models (GAMs). LOI, magnetic susceptibility, and XRF elements Fe and Ca, and Si/Al and Si/Ti ratios were used as predictor variables. Fe and Ca were identified as key elemental variables by performing a principal component analysis on the XRF data (Supporting Information Fig. S16). Si/Al and Si/Ti ratios were used as proxies for biogenic silica (Supporting Information Figs. S17-S19) (Meckler et al. 2013;Balascio et al. 2017). Due to the difference in sampling resolution, geochemical values were calculated for the chironomid sample dates using the predict function in R (Supporting Information Figs. S20 and S21). The relationship between mean July temperature and the chironomid ecological metrics were estimated using GAMs for the period 1901-2018. Local historical (mean July) temperature observations were available from 1901 to 2018 from the CRU TS4 data set to the nearest latitude-longitude half degree grid (Harris et al. 2020). Temperature values were calculated for the chironomid sample dates using the predict function in R (Supporting Information Fig. S22). GAMs were run using the gam function from the mgcv package in R (Wood 2011;Wood et al. 2017). For full model outputs, see Supporting Information Tables S3-S20.

Simulated changes in taxon connectivity
The taxon-richness metric clearly distinguished between taxon loss, gain, and replacement (Fig. 3, left-hand columns). However, it had little sensitivity to structural change as simulated by changes in taxonomic connectivity. Taxon richness did not show unique signatures for non-random change in connectivity to distinguish from random change.
Under taxonomic loss, beta diversity demonstrated more sensitivity to structural change than taxon richness or skewness (Fig. 3, scenarios 1-3). Beta diversity increased in all three scenarios indicating a rise in dissimilarity with taxon loss. Beta diversity became more variable with the loss of weakly connected taxa as the networks decreased in size (scenario 1). The loss of strongly connected taxa produced a gradual increase in beta diversity (scenario 2), distinct from the loss of randomly connected taxa. The random loss of taxa, irrespective of connectivity, produced the greatest increase and variation in beta diversity (scenario 3).
Under taxonomic gain, skewness had more sensitivity to structural change than taxon richness or beta diversity (Fig. 3, scenarios 4-6). The gain of weakly connected taxa had a unique signal, with a small, abrupt decrease in skewness after which skewness remained relatively unchanging (scenario 4). The gain of strongly connected and randomly connected taxa produced large decreases in skewness. There was a small rise in skewness as strongly connected taxa continued to join the network (scenario 5). Skewness continued to decline gradually as randomly connected taxa continued to join the network (scenario 6).
Under taxonomic replacement, both beta diversity and skewness were required to see the unique signatures of the type of structural change (Fig. 3, scenarios 7-11). Beta diversity increased and skewness experienced a small decrease as weakly connected taxa joined the network structure (scenarios 7 and 10). As existing taxa were replaced with strongly connected taxa, beta diversity was highly variable, with little directional trend, and skewness strongly decreased (scenarios 8 and 9). In the random replacement scenario, beta diversity was most variable at the ends of the model run, and skewness decreased with taxon replacement (scenario 11), most similar to the replacement of weakly connected taxa with strongly connected taxa (scenario 9).
Overall, the simulation models indicated taxon richness was not sensitive to changes in taxon connectivity, and beta diversity was more sensitive to changes in connectivity in the taxon loss and replacement scenarios. Skewness was more sensitive to changes in connectivity associated with taxonomic P a r a t a n y t a r s u s p e n ic il la t u s -t y p e P a r a t a n y t a r s u s u n d iff T a n y t a r s u s lu g e n s -t y p e T a n y t a r s u s n o s p u r T a n y t a r s u s p a ll id ic o r n is -t y p e T a n y t a r s u s lu

P s e c t r o c l a d i u s s o r d i d e l l u s -t y p e P a r a k i e ff e r i e l l a t r i q u e t r a -t y p e C h i r o n o m u s a n t h r a c i n u s -t y p e P a r a t e n d i p e s a l b i m a n u s -t y p e
T a n y t a r s u s g l a b r e s c e n s -t y p e T a n y t a r s i n i u n d i f T a n y t a r s u s l u g e n s -t y p e T a n y t a r s u s n o s p u r P e n t a n e u r i n i u n d i f P a g a s t i e l l a C l a d o p e l m a P a r a t a n y t a r s u s u n d i ff C r i c o t o p u s -t y p e P

P s e c t r o c l a d i u s s o r d i d e l l u s -t y p e C h i r o n o m u s a n t h r a c i n u s -t y p e P a r a c h i r o n o m u s v a r u s -t y p e
T a n y t a r s u s g l a b r e s c e n s -t y p e T a n y t a r s i n i u n d i f T a n y t a r s u s l u g e n s -t y p e T a n y t a r s u s n o s p u r

Assemblage and structural change in chironomid communities
The chironomid assemblage at Round Tangle lake (Fig. 4) was dominated by Sergentia coracina-type, Paratanytarsus undifferentiated (no mandibles), Tanytarsus no spur, Tanytarsini undifferentiated, Orthocladius consobrinus-type, Paracladius, and Zalutschia type B. The cluster analysis indicated that the greatest change in taxonomic composition occurred between the most recent sample and the rest of the record. The most recent sample (Zone 2) indicated an increase of S. coracina-type, Corynocera oliveri-type, Cricotopus cylindraceus, Diplocladius, Limnophyes-Paralimnophyes, and Thienemannimyia, and a decrease in Paratanytarsus undifferentiated, Cricotopus intersectus, O. consobrinustype, Paracladius, and Procladius.
All three data sets indicated variation in ecological space across the DCA axis 1 (Fig. 7, upper panels). Changes in the assemblage composition zones, as indicated by dendrographs (Figs. 4-6), plotted separately across the DCA axis 1 (light toned symbols indicate younger samples and dark toned symbols indicate the older samples). CCA plots indicate that changes in the organic (LOI) and minerogenic content were likely to have influenced the chironomid assemblages (Fig. 7, lower panels). The historical temperature data corresponded with LOI in a CCA applied to the data covering the period c. 1901-2018, suggesting the change in productivity corresponded with changing temperature (Supporting Information Fig. S23).
The three lakes had unique faunal assemblages, as expected for lakes in different locations with different water conditions (see Supporting Information Table S1), and thus had different changes in taxa present in response to the changing climate/ environmental conditions. In this paper, we are interested in the relative changes in community composition and structure, as represented by rarefaction, beta diversity and skewness, and whether these changes were comparable despite the differing faunal assemblages.
Rarefaction increased in the upper assemblage samples in all three lake records; however, there was no detectable change in the median rarefaction values across the assemblage zones in any of the records (Fig. 8). GAMs indicated relationships between the sedimentological variables and rarefaction in Round Tangle and East Cobb lakes. In Round Tangle, Fe and Ca explained 47.1% of the deviance in rarefaction (p value < 0.005). In East Cobb, LOI explained 30.9% of the deviance in rarefaction (p value < 0.005).
Beta diversity increased detectably in East Cobb, indicating a rise in turnover in association with the identified assemblage zones (Fig. 8). Beta diversity did not change discernibly across the assemblage zones in Round Tangle or Little Lost. The GAMs indicated a strong relationship between LOI and beta diversity, with LOI explaining 90.3% of the deviance in beta diversity in Round Tangle and 63.2% of the deviance in East Cobb (p values < 0.001). GAMs also indicated a strong relationship between the Si/Al and Si/Ti ratios and beta diversity, explaining 91.2% of the deviance in beta diversity in Round Tangle (p values < 0.001 and 0.01) and 57.5% of the deviance in East Cobb (p values 0.01 and 0.02). A GAM indicated a strong relationship between July temperature and beta diversity in Round Tangle, with July temperature explaining 97.1% of the deviance in beta diversity for the period 1901-2018.
Round Tangle had the least negative skewness values, with East Cobb and Little Lost having more negative skewness values ( Fig. 8; Supporting Information Fig. S24). There was little directional change in the skewness values for Round Tangle and East Cobb, although in East Cobb, the Norwegian skewness model suggested a relative decline in the latter half of the record (Fig. 8). In Little Lost, skewness decreased in the recent samples, for both skewness methods. Skewness trends indicated different fluctuations in all three lakes based on the filtering technique ( Fig. 8; Supporting Information Figs. S14, S24). Generally, skewness values calculated using the Norwegian calibration data set were less negative than skewness values calculated using the North American calibration data set (Supporting Information  Fig. S24). The three records showed no detectable change in assemblage structure in association with the assemblage zones identified by the dendrographs.

Discussion
High-latitude lakes have been exposed to a wide range of stressors over recent centuries , driving changes in freshwater diatom and Pediastrum algae communities (Smol and Douglas 2007b;Woelders et al. 2018). Here, we demonstrate that recent chironomid communities in Subarctic Alaska have also experienced compositional change in relation to climatic and environmental stress. Comparison between the metric trends in the simulated and empirical data sets suggests that the empirical data sets have maintained some degree of structural stability under the current levels of climate and environmental change.
Network structure as an indicator of ecosystem stability Network structure can indicate the level of stress a community might be experiencing, through the proportion of strongly and weakly connected taxa (Wang et al. 2019). The network models in this paper suggest that taxa joining a network have a greater influence on the network structure than taxa leaving a network. As taxa join the networks, skewness declines; with the arrival of weakly connected taxa causing a smaller decrease in skewness (scenarios 4, 7, and 10) and the arrival of strongly connected taxa causing a larger decrease in skewness (scenarios 5, 8, and 9). These trends are evident in the gain and replacement scenarios, suggesting that the arrival of new taxa has greater influence on the network structure than the loss of taxa. The loss scenarios support this, indicating that the loss of taxa simulated gradual changes in skewness, irrespective of the connectivity of the lost taxa (scenarios 1-3). In the hypothesized models in Mayfield et al. (2020) where connectivity was not accounted for, the arrival of taxa also simulated a decrease in skewness values in scenarios where taxon richness increased linearly and where taxon richness increased overall with turnover in taxon type. Thus, the models used in this study indicate that the relative connectivity of the taxa joining a network is key to understanding the effect of these additional taxa on the network structure. It is likely that the connectivity of a joining taxon is determined by the heterogeneity/homogeneity of the environment the taxon is joining. Under reduced exogenous forcing, systems dominated by weakly connected taxa are likely to have relatively high diversity and turnover rates, where individual taxa react to changes in environmental conditions (leave or join the ecosystem) at their own rate (van Nes and Scheffer 2005;Scheffer et al. 2012). This is unlikely to have a strong effect on the network structure, as suggested by modeled scenario 7. However, increased environmental stress from exogenous forcing has been associated with the increased prevalence of strongly connected taxa in analyses of community nestedness (Doncaster et al. 2016). Similarly, lake environmental homogenization, due to eutrophication, corresponded with abrupt declines in the network skewness of diatom communities (Wang et al. 2019). Increased stress can cause environmental homogenization and a loss of microhabitats increasing the number of interactions between existing taxa, and perhaps increasing the number of interactions a newly arrived taxon is likely to construct. Understanding community structure, through analyses such as network skewness, may provide earlier detection of ecosystem instability, increasing the window of opportunity to manage or reduce shifts in ecological state (Hughes et al. 2013).

Relative ecosystem stability in changing Subarctic Alaskan lakes
Changes in the taxonomic composition and sedimentology in the lakes presented here suggest that these lakes have experienced recent climatic or other environmental change. A decrease in cold, oligotrophic taxa, such as Co. ambigua and Co. oliveri-type, and rise in temperate-preferring taxa, such as Ch. anthracinus-type, Mi. pedellus-type, and Po. nubeculosum-type, in East Cobb lake indicated a direct response to temperature change (Brooks et al. 2007). In Round Tangle, there was a rise in the presence of littoral taxa, such as Diplocladius, Limnophyes-Paralimnophyes, and Thienemannimyia-type, which could indicate lake shallowing (Cranston 1983;Fittkau and Roback 1983;Massaferro and Brooks 2002;Brooks et al. 2007). Lake shallowing and shrinking have been recorded across Subarctic and boreal Alaska, attributed to climate-driven processes, such as increased drainage from permafrost melting and increased evapotranspiration from warmer and longer summers (Riordan et al. 2006;Jepsen et al. 2013). Thienemannimyia-type and Limnophyes can also be found in streams (Fittkau and Roback 1983;Brooks et al. 2007), thus, a rise in these taxa could indicate increased catchment activity or inwash. However, a decrease in concentration of the geochemical composition within the recent sediments at Round Tangle (Supporting Information Fig. S9) suggests increased inwash was not the primary reason for the increased presence of littoral taxa (Schlolaut et al. 2014;Davies et al. 2015;Plaza-Morlote et al. 2017).
The trends in rarefaction and beta diversity are congruent with changes in climate and/or environmental conditions. Increased chironomid richness can be associated with warming climates; cold environments are often taxon poor (Brodersen and Anderson 2002;Brooks and Birks 2004) as temperature directly affects chironomid physiological and biochemical processes such as respiration, development, and growth (Eggermont and Heiri 2012). All three lakes indicated an overall rise in the trend of rarefaction over the last c. 100 yr. The changes in beta diversity in the chironomid communities were strongly related to changes in biological productivity, likely linked to warmer temperatures (as identified by the CCA and GAM outputs). In Little Lost, there was a rise in the abundance of Di. nervosus-type and Ablabesmyia, which are often associated with macrophytes (Brodersen et al. 2001;Langdon et al. 2010;Watson et al. 2010), and a decrease in Pagastiella, which can be sensitive to increases in eutrophication (Moller Pillot and Buskens 1990). In East Cobb, there was an increase in the abundance of Di. nervosus-type (associated with macrophytes) and Po. nubeculosum-type, which often occurs in eutrophic lakes (Brooks et al. 2007). Increased biological productivity in lakes has been linked with warming in the Arctic (Michelutti et al. 2005). Thus, it is likely that the assemblages in this study, relating to increased productivity and potentially lake shallowing, are showing signs of ecosystem change related to climate change. The relative increase in rarefaction and turnover within the recent chironomid samples suggest that the empirical records are more comparable to the hypothetical gain or replacement scenarios.
Skewness was predominately negative across all three lakes, indicating relatively high taxon connectivity. This suggests some degree of community stress; Mayfield et al. (2020) found that the majority of high latitude chironomid communities were experiencing structural stress, likely due to their extreme environments. There was little discernible change in community structure across the records; however, in the Little Lost record, both skewness methods indicated more negative skewness in the upper samples suggesting a higher proportion of strongly connected taxa. More negative skewness is associated with increases in taxon richness and beta diversity, congruent with modeled scenario 5, where strongly connected taxa join the existing network. Thus, the Little Lost record could be showing signs of structural instability; however, the Little Lost record is short (10 samples long), making it difficult to distinguish the trend from background variation.

Possible future trajectories and resilience of Subarctic chironomid communities
The chironomid communities presented in this paper show signs of climate-and environment-related community change. However, the communities have retained a degree of structural stability, as indicated by the lack of detectable change in skewness. Cold stenothermic taxa are adapted to thrive in harsh environments, perhaps providing Arctic and Subarctic ecosystems with extra resilience to changing conditions (Stewart et al. 2013). Chironomids can have relatively large trait variability or plasticity (Serra et al. 2017) and opportunistic feeding behaviors (Lee et al. 2018), further promoting a degree of resilience within chironomid communities. Theoretically, highly connected networks can provide a temporary resilience to stress. However, they have less capacity to respond gradually to increased stress as individual taxa cannot respond to stress in isolation but rather in connection with other taxa within the ecosystem (van Nes and Scheffer 2005;Walther 2010). This increases the likelihood of systemic critical transitions; a sudden and system-wide shift in taxa (Scheffer et al. 2012). Thus, these highly connected taxa may be vulnerable to abrupt and system-wide structural change under continued climate and environmental stress. An understanding of the structural stability of aquatic systems, through methods such as skewness, can help anticipate ecosystem structural change and thus enable the implementation of management techniques to reduce the likelihood of system-wide structural changes.
The rise in rarefaction and turnover in the chironomid communities was associated with a rise in biological productivity (LOI %, Si/Al and Si/Ti ratios) and possibly lake shallowing. The rise in productivity could have increased the availability of microhabitats within the lake, for example, through an increase in macrophyte presence (Langdon et al. 2010). Aquatic vegetation was present in all three lakes on sampling. This rise in microhabitat diversity, or heterogeneity, may explain the lack of definitive change in network skewness. Hayden et al. (2017) also found a peak in macroinvertebrate community in mesotrophic lakes in the Subarctic. Mayfield et al. (2020) found that initial warming may provide some ecological release for high latitude lakes with extreme environments with taxon-poor cold-adapted assemblages. Warmer temperatures can increase biological productivity, increasing lake habitat diversity and, thus, decreasing taxonomic interactions and structural stress. It is possible that the Round Tangle and East Cobb chironomid assemblages are in this initial warming stage, with increased richness and turnover, and relative structural stability provided by the increased lake biological productivity and macrophyte presence. However, continued temperature stress, and potentially the addition of secondary drivers of stress, may cause ecological degradation through habitat homogenization and increase structural stress (Mayfield et al. 2020). Little Lost lake has experienced higher temperatures than Round Tangle and East Cobb lakes (Fig. 2), thus it may have experienced greater temperature stress, perhaps explaining why Little Lost demonstrates a clearer negative trend in skewness. Continued monitoring of Little Lost lake may provide supporting evidence of the declining structural stability of this lake, while the application of these metrics to a greater number of lakes in similar environments may provide greater indication as to whether these metrics can act as an early warning signal of structural instability.
As the climate continues to warm (Smith et al. 2015) with associated increases in eutrophication (Douglas and Smol 1999;Antoniades et al. 2011;Hessen et al. 2017), the amount of stress high latitude lakes are exposed to is likely to increase. It is thought that high rates and magnitudes of climate change are required to alter aquatic community structures (Mayfield et al. 2021). Thus the lack of structural change in the chironomid communities presented here suggests that the current level of climate and environmental change experienced by the lakes has not reached a high enough rate or magnitude to drive structural change. Alternatively, it is possible that the current chironomid communities are demonstrating a degree of relative resilience; Scheffer et al. (2012) suggest that highly connected networks have a temporary resilience to stress. The lack of trend in the skewness values in the lakes presented here may indicate some resilience to the amount of stress currently experienced by these lakes. With current predictions indicating that high latitude regions are likely to experience more warming (Overpeck et al. 1997;Post et al. 2009;Pithan and Mauritsen 2014), it is likely that high latitude lakes will experience greater levels of stress, causing further taxonomic turnover as cold adapted taxa are lost and temperate or more (meso)eutrophic taxa arrive. Decreases in skewness have already been recognized in contemporary lake assemblages driven by eutrophication and ecosystem homogenization (Wang et al. 2019) and in subarctic lakes exposed to warmer temperatures and possible secondary variables (Mayfield et al. 2020). Thus, the increasing trends in rarefaction and turnover in Round Tangle and East Cobb lakes, and increasingly negative skewness values in Little Lost lake could be suggestive of further compositional change under continued climate and environmental stress.

Taxonomic resolution and network skewness
Skewness produced different trends, based on usage of the North American and Norwegian calibration data sets, with little correlation. This suggests that the filtering of the fossil data sets can affect skewness, most likely because the grouping process obscures the interactions between taxa. In the North American calibration data set, Cricotopus-Orthocladius speciestypes are grouped and the majority of Tanytarsini taxa are combined in a sub-tribe Tanytarsina-other group, including rare/ common taxa and warm/ cold taxa. For example, Tanytarsus lugens is a cold-preferring taxon and commonly occurs with other cold-water taxa, while Tanytarsini glabrescens is more commonly a warm-indicator taxon (Watson et al. 2010;Brooks et al. 2012). It is unlikely for these Tanytarsus species to occur simultaneously, however, they are grouped in the North American data set, changing the co-occurrences of these taxa with others and likely effecting skewness. In the Little Lost and East Cobb records, changes in skewness were more pronounced when calculated using the Norwegian calibration data set. The Norwegian calibration data set has a higher taxonomic resolution, and thus, perhaps provides a more accurate representation of structural change.

Conclusions
This study demonstrates the current chironomid community structure of three Subarctic lakes in Alaska using taxon richness (rarefaction), beta diversity and network skewness. The sensitivity of these metrics to changing taxonomic connectivity was tested using model-simulated taxonomic networks. These simulations indicated that taxon richness was not sensitive to connectivity, beta diversity was more sensitive to taxon loss, and skewness was more sensitive to taxon gain. Beta diversity and skewness were both required to understand changes in network structure under taxon replacement scenarios. The taxon gain and replacement scenarios indicated that the arrival of strongly connected taxa simulated a greater decrease in skewness than the arrival of weakly connected taxa. In the empirical records, increases in the rarefaction and beta diversity trends and changes in the taxonomic composition indicated that Round Tangle and East Cobb lakes have experienced recent climatic or other environmental change, without discernible changes in community structure, suggesting some degree of structural resilience to the current levels of climatic and environmental stress. In Little Lost lake, decreases in the skewness trend were indicative of structural instability, a possible sign of ecological damage.
Our definition of ecosystem resilience was the capacity of an ecosystem to absorb disturbance and retain the same ecological function, without shifting into an alternative state. In these lakes, we find that in the main structural stability, and hence resilience is maintained for now at least, although in Little Lost Lake, we find indicators of structural instability, and hence a potential loss of resilience to climatic and environmental stress. Applying these ecological metrics to a greater number of lakes in similar environmental conditions may provide greater indication as to whether these metrics can act as an early warning signal of structural instability.