Does glacial retreat impact benthic chironomid communities? A case study from Rocky Mountain National Park, Colorado

The aim of this study was to determine which environmental variables are responsible for modern benthic chironomid distributions in a glacial setting. The chironomid communities from nine alpine lakes were assessed, and forty-three individual taxa were extracted and identified. Surface water temperature and nitrate were strongly and negatively correlated (−0.82, p = 0.007), suggesting that glacial meltwater (the driver that explains both surface water temperature (SWT) (°C) and nitrate (NO3 + NO2-N)) is the environmental variable that explains the most variance (15%). On average, lakes receiving glacial meltwater were 2.62 °C colder and contained 66% more NO3 + NO2-N than lakes only receiving meltwater from snow. The presence of taxa from the tribe Diamesinae indicates very cold input from running water, and these taxa may be used as a qualitative indicator species for the existence of glacial meltwater within a lake catchment. Heterotrissocladius, Diamesa spp., and Pseudodiamesa were present in the coldest lakes. Chironomus, Diplocladius, and Protanypus were assemblages found in cold lakes affiliated with the littoral zone or alpine streams. The modern benthic chironomid communities collected from the alpine of subalpine lakes of Rocky Mountain National Park, Colorado, represent a range of climatic and trophic influences and capture the transition from cold oligotrophic lakes to warmer and eutrophic conditions.


Introduction
Globally, the remote lakes found at high elevations act as sentinels of change and are among the first bodies of water impacted by climate change [24]. As anthropogenic warming continues, the physical and chemical limnology of these lakes will change. As a result, the structure, function, and composition of aquatic ecosystems in subalpine and alpine settings will be significantly transformed [71]. Recent studies suggest that the range of invasive fish species and invertebrates will expand due to their ability to move to higher elevation lakes via warmer creeks with stable environments [34]. This movement will impact natural communities that exist in alpine lakes by altering food web and predator/prey dynamics [34]. Increases in terrestrial vegetation related to the upslope movement of timberlines will supply higher amounts of dissolved organic carbon (DOC) to subalpine and alpine lakes, which will increase productivity and potentially affect the diversity of high lake ecosystems [44].
Many of the alpine lakes in the western USA receive cold meltwater emanating from small cirque and rock glaciers. These alpine environments are sensitive to regional climatic change [16,43]. Hall and Fagre [22] modeled glacial retreat for Glacier National Park and concluded that glaciers would disappear from the landscape by 2030. Recent work suggests that local topographic effects may buffer against regional warming and glacial extinction may be delayed by 50-250 years [5,13,20,64]. Studies indicate that glacial meltwater affects alpine hydrology, chemistry, and the turbidity of high alpine lake water [74]. Temperatures in the western USA have steadily increased over the past few decades and have amplified the rate at which glaciers and permafrost are melting in alpine areas [14]. The addition of this cold, silt-enriched water into alpine lakes will impact the timing of lake stratification [29], increase turbidity, decrease optical transparency [53], and alter the water chemistry [41] of these alpine lakes. Studies indicate that lakes receiving glacial meltwater have up to 200× more nitrogen than those lakes that only receive snowmelt [9,47,50,74]. Atmospheric nitrogen has been accumulating on the surface of glaciers for decades. This nitrate is then added to lake system with the onset of glacial retreat. The additional input of nitrogen to N-limited lakes promotes enhanced primary productivity, such as algal blooms [21,73]. While many studies have endeavored to elucidate the physical and geochemical changes within these systems, little work has addressed how biotic communities in alpine lakes will respond to the increased flux of glacial meltwater in the short term [74]. Limited studies have documented that diatom communities in many alpine lakes are shifting from large filamentous diatoms (i.e., Aulacoseira taxa) to Cyclotella spp. as a result of longer growing seasons, increased stratification, and decreased ice cover [12]. In addition, Fegel et al. [16] found changes in microbial communities due to the geochemical changes present in the headwaters emanating from rock glaciers and glaciers from the Cascade Mountains, the Sierra Nevada, and the Rocky Mountains of the western USA.
Chironomids, or midges, are some of the most abundant insects found in freshwater ecosystems and are useful to study changes in temperature, pollution, and dynamic system changes [8,17,60,79]. Midges occupy several trophic levels in aquatic ecosystems and therefore play a vital ecological role in lakes [60,79]. However, the ecological understanding of chironomid distribution is poorly understood [15].
A chironomid life cycle includes several stages that begins as an egg mass deposited on the surface of the water by an adult chironomid. As the eggs hatch, chironomids erupt in their first larval state and mainly persist as benthos on the floor of the lake. In this state, the chironomid has a maggot-like form and a chitinous head capsule that is shed three more times as individuals grow [60,79]. Eventually, the larval chironomid reaches the pupae stage and rises through the water column of the lake. This stage is abrupt and leads to metamorphosis from pupae to an adult fly that emerges from the lake [60,79]. Thus, the survival of chironomid egg masses is largely influenced by surface water temperature [72], the larval stage is influenced by bottom water temperature, and adult flies exist in environments dominated by air temperature [15]. Chironomids have been used as a biological proxy to model both surface water temperature [7,61,83] and air temperature [23,26,35,38] based on the assumption that a strong relationship exists between surface water temperature and air temperature. Eggermont and Heiri [15] caution that multiple factors such as depth, thermal stratification, and glacial melt may impact the relationship between air and surface water temperatures. Understanding the dichotomy between air and water temperatures in chironomid ecology is imperative for future chironomidbased paleoclimate studies.
To date, the studies that have assessed the response of midges to glacial melt in alpine settings focus on montane streams [36,37,46,69]. These studies indicate that cold water obligate chironomid communities are found in glacial meltwaters; however, there remains a paucity of studies documenting the response of midges to glacial melt in lacustrine settings. Information extracted from lake sediment can be used to develop baseline limnological information against which future changes can be compared [76]. This research is vital due to the very narrow window of time that is left for studies that examine glacial retreat due to the projected demise of alpine glaciers, especially those present in the western USA [73].

Study area
The Colorado Rocky Mountains possess the most southern-reaching alpine glaciers currently still active in the USA. While many areas present in the Northern and Central Rocky Mountains have shown pronounced ablation rates for alpine glaciers [2,49], the behavior of glaciers in the Southern Rocky Mountains is quite different and the rate of glacial retreat is much slower in this region relative to regions to the north, such as Glacier National Park [27,55,65,66], and makes the Front Range a critical location for monitoring glacier change [27].
Rocky Mountain National Park (RMNP) is in the northern portion of the Front Range and is home to 30 glaciers [27]. The glaciers in RMNP straddle an elevational range that includes regional timberline (3500 m Above Sea Level or ASL) and lies between 3416 and 4068 m ASL. Most glaciers are found on the eastern side of the Continental Divide and occupy north-to east-facing cirques [39]. Snow accumulation is frequently redistributed into these cirques by strong westerly winds and avalanching [27,54,86]. The local topographic shading evident on the eastern side of the Continental Divide also has strong control over ablation rates and may account for the highly irregular ablation-altitude gradients evident in the Front Range [27,54].
Comparing the chironomid assemblages found in fives lakes located in glacial basins to four lakes in snow basins (with variations in elevation, geology, and vegetation controlled for) enabled an assessment of the relative role meltwater plays in shaping chironomid communities found in the alpine lakes located in RMNP. Lakes associated with glacial ice bodies were selected based on Fountain et al. [18] and the National Geographic Rocky Mountain National Park topographic map (#200). All lakes were found in areas that consisted of igneous Proterozoic diorites and granites that intrude into ancestral metamorphic Proterozoic biotite gneisses, migmatites, and schists [33].
The collection of ten short lacustrine sediment cores occurred during the late summers of 2015 and 2016 (Fig. 1). Study sites were chosen to be at the approximate elevation with similar vegetation and geology. Cony (CNY) and Pipit (PIP) lakes were the highest elevations and were sampled at 3509 and 3479 m ASL, respectively (see Table 1). These lakes are located in rocky cirques above timberline adjacent to the continental divide. Cony Lake receives glacial melt from an unnamed glacier, whereas Pipit only receives melt from annual snowfall. Hutcheson Lake (HCH) (3413 m ASL) receives glacial input from Cony Lake, which lies immediately above the catchment. Grasses that are typical of alpine-tundra surround the lake. The snowmelt-fed Falcon Lake (FAL) lies at 3371 m ASL and is located in a small rocky cirque with small patches of krummholz. Box (BOX) and Eagle (EGL) Lakes lie at 3274 and 3298 m ASL, respectively, and are located at timberline. Eagle Lake receives glacial meltwater from Moomaw Glacier, whereas Box Lake only receives annual snowmelt. The remaining lakes all lie below timberline and are in subalpine forest predominately composed of Pseudotsuga menziesii (Douglas fir) and Picea engelmannii (Engelmann spruce). Black (BLK) Lake, located at 3237 m ASL, is the deepest lake sampled at 21.2 m and receives glacial melt. BLK also has steep scree slopes SSW and S of the lake. Thunder (THD) Lake lies at 3225 m ASL and is on a much gentler slope than BLK. The lake is surrounded by forest almost to the edge of the lake and only receives annual snowmelt. Odessa (ODS) and Spruce (SPR) Lakes were the lowest in elevation of the sampled lakes. ODS (3051 m ASL) receives glacial input from three unnamed glaciers located higher in the catchment. Spruce Lake (2947 m ASL) was unlike any other lake sampled. It was only 1 m deep and had tall grasses throughout the entire bed of the lake. It also contained no chironomid subfossil remains and thus was eliminated from the analysis.

Field methods
The sediment cores were collected from the approximate center of the study lakes using a gravity corer, deployed from a small, two-person inflatable raft that allows recovery of lake surface sediment with minimum disturbance of the mud-water interface. The cores were typically 20 cm long and represented approximately 150 years of deposition. Observations regarding the stratigraphy and color of each core were recorded in a field notebook. Each core was sectioned into 0.25-cm intervals and placed into Whirl-paks®. A Yellow Springs Instrument (YSI®) Professional Plus was used to collect a suite of limnological variables, such as temperature, pH, and specific conductivity. Water samples were collected from the center of each lake and submitted to the Center for Applied Isotope Studies (CAIS) at the University of Georgia for analysis of analytes for nitrogen, phosphorus, chlorophyll α (chl a), as well as nutrients for chloride and sulfate. The cores were transported to the Environmental Change Lab at the University of Georgia in coolers after the end of each field season. Total carbon (%) of dry bulk sediment was analyzed using EA-IRMS at CAIS. Distance (m) is the measured distance from the lake to the terminus of the glacier in the lake catchment. The glacial index (GI) is a measure of environmental harshness and is an index of glacial influence following Jacobsen and Dangles [32]. The GI was calculated as GI = √ area distance+ √ area (for area > 0). The glacial coverage in the catchment (GCC) was also calculated [31]. The area of each catchment as well as those for each glacier was determined using the GLIMS Glacier Database in ArcGIS [67]. Forty environmental variables were collected in total ( Table 2). Physical variables included elevation, lake depth, Secchi disk depth, mean July air temperature, surface water temperature, bottom water temperature, the GI, GCC, and the distance from the glacier. Geochemistry variables that were sampled include pH, specific conductivity, dissolved oxygen, dissolved organic carbon, dissolved inorganic carbon, total phosphorus as PO 4 -P, and NO 3 -N + NO 2 -N among others (Table 2). However, an additional 26 elements collected from lake water were below detection and removed from the analysis.

Laboratory
Chironomid extraction procedures followed the protocol established by Walker [84]. Bulk sediment samples were soaked in an 8% KOH solution and heated to 40 °C for a minimum of 30 min. The solution was then sieved through a 95-μm grade mesh screen using distilled water to eliminate any remaining KOH residue. The material remaining on the screen was transferred into a beaker with distilled water. The resulting residue was then poured into a Bogorov counting tray and sorted using a stereoscope at 40X. The sub-fossil chironomid head capsules extracted from the residue were permanently mounted on glass slides using Entellan®. This process was repeated until a minimum of 50 head capsules was recovered from each sample following the advice of Heiri and Lotter [25]. A Nikon Eclipse E100 (×100) microscope was used for taxonomic determination of the midge remains. The taxonomic keys by Brooks et al. [8] and Andersen et al. [1] were instrumental in the identification of midge taxa. A total of 542.5 head capsules were collected and counted from the top 0.50 cm of sediment of each lake to assess the modern distribution of chironomid communities (mean = 60.28, maximum = 88, minimum = 43.5). Fortythree taxa were identified from the modern sediment.  However, if less than 2% of a particular taxon was represented and they were present in fewer than 2 lakes, they were removed from statistical analyses following Quinlan and Smol [62]. Thus 30 taxa were used in analysis. Spruce Lake contained zero chironomid remains and was removed from the analysis.

Data analysis
Detrended canonical analysis (DCA) is an indirect ordination technique that is useful in the exploration of taxa data collected from the lacustrine sediment. Chironomid taxa possess the highest abundances in environments that maximize their preferred habitats. Abundances begin to decline or disappear as they become farther removed from their preferred environment. Due to these characteristics, ecological data typically possess "a modal relationship to their ecological gradients" [28]. DCA assumes that the data have a unimodal distribution. The chironomid assemblage data followed a Poisson distribution and satisfied this assumption. The data were square-root transformed to shorten the distribution and to make the data homoscedastic. The effect of rare taxa was down-weighted to dampen their effects on the ordination. DCA is used to determine whether a linear, e.g., redundancy analysis (RDA), or unimodal, e.g., canonical correspondence analysis (CCA), model should be used to understand which environmental variables explain the most variance in the distribution of chironomid communities [77]. The length of the first DCA axis was 2.48, and RDA was chosen as the appropriate model to assess the relationship between chironomids and environmental variables. Redundancy analysis (RDA) was used to assess which environmental variable explains the most variance in species distribution [78,88].
RDA is used to extract and summarize variation in the chironomid taxa data that can be explained by environmental variables [88]. RDA models were developed using forward selection combined with Monte Carlo permutation tests (p < 0.05, 999 permutations) in order to identify the environmental variables that most likely explained the distribution of modern chironomid communities. This method, combined with permutations, was used to determine the statistical significance (p < 0.05) of each environmental variable, as well as demonstrating the amount of variance that each variable accounted for [6,77]. Forty environmental variables were collected. Physical variables included elevation, lake depth. However, 26 elements collected from lake water were below detection and removed from the analysis. Variables that covaried, such as dissolved oxygen (%) and dissolved oxygen (mg/L), were examined, and only one representative variable was used. Twenty-five remaining environmental variables were assessed for linearity, and specific conductivity, lake depth, and Secchi disk depth were log-transformed to ensure that homoscedasticity assumptions were met ( Table 1). All statistical analyses were performed using the open-source platform R (version 3.6.1) (R Development Core Team [63], http:// www.R-proje ct. org). DCA and RDA were implemented in the vegan package [51].

Water chemistry
In general, the lakes sampled were relatively deep, and all were over 5 m with the exception of Hutcheson Lake (3.1 m). Lake depths ranged from 3.1 to 21.2 m deep, with an average of 10.42 m. The surface temperature for lake water was variable and ranged from 8.2 to 13.9 °C. This wide range in surface water temperature (SWT) is of note as these lakes are not found on a particularly long elevational gradient (458 m). All lakes were open basins and received input from alpine streams as well as having outlet streams. The temperature profiles for the shallowest lakes that received glacial meltwater (i.e., Odessa and Hutcheson Lakes) showed no sign of thermal stratification and only varied by ≤ 1.4 °C. Pipit and Falcon Lakes were deeper (10.4 m and 8 m, respectively) but only had temperature profiles that varied by 1 °C from the surface water to the bottom of the lake. These lakes are fed only by annual snowmelt. Cony, Eagle, Box, Black, and Thunder Lakes all possessed temperature profiles indicative of thermal stratification with an epilimnion thickness of 6-8 m.
The trophic class of each lake was evaluated following Carlson's Trophic State Index [11]. The Carlson Index uses three independent variables of aquatic biomass that includes Secchi depth (SD), total phosphorus (TP) from the epilimnion, and chlorophyll α (Chl) ( Table 2). However, Horne and Goldman [30] warn that trophic classifications are idealized concepts and that real-world examples are more varied. For this study, lakes were classified if two of the three variables fell within range of a specific trophic level. The results indicate that only Pipit Lake may be considered oligotrophic (Tables 1 and 2). Total P was below  The only eutrophic lake sampled was Thunder Lake. TP was high at 52.08 ppb, and the SD was 2.2 m. The remaining seven lakes are classed as mesotrophic, which are typically lakes with an intermediate level of productivity [30].

Chironomids
Chironomus spp. was the dominant taxa present and comprised 24.3% of the total chironomids recovered from all lakes (Fig. 2). Chironomus is eurythermic and is known as a "blood worm" as it emits a red color due to the hemoglobin it produces [56]. It is mostly found in the profundal zone (i.e., the deepest zone) of lakes and can tolerate low levels of oxygen and or short periods of anoxia for this reason [8,85]. It is opportunistic and is often found in lakes undergoing environmental change as it is an early colonizer [8]. Corynocera oliveri type was the second most abundant taxa (10.9%). This taxon is typically found in the muddy substrate of cold lakes [1,8]. Porinchu and Cwynar [57] documented the presence of these insects with regard to timberline in Siberia. They found that C. oliveri was found typically in the colder lakes located above timberline. While this is true of the assemblages collected from Pipit, Hutcheson, Eagle, and Box Lakes, C. oliveri had higher relative abundances from lakes below timberline (Thunder Lake and Odessa Lake). Heterotrissocladius spp. (6.9%) is a very common taxon in all lakes collected and is typically found in the profundal of cold oligotrophic lakes that are well-oxygenated [1,8,80]. Procladius (3.7%) was also one of the most prevalent taxa. However, the relative abundance of Procladius is very high in Pipit Lake and much lower in every other lake. Procladius is very common in lakes that are classified as mesotrophic and eutrophic and is typically associated with the profundal zone [8]. Sergentia (3.52%) is typically found in relatively deep, mesotrophic lakes. This taxon is found in all sampled lakes except for Cony and Eagle Lakes. The presence of both Chironomus and Sergentia indicates early colonization is occurring, and a transition from an oligotrophic state to a mesotrophic trophic state is in progress. The presence of taxa from the tribe Diamesinae (including Diamesa spp., Pseudodiamesa, and Protanypus) is of particular interest in that the remains of these taxa are extremely rare in lake sediment and poorly studied [8,56,82]. Recent studies of chironomids in alpine streams find that the presence of Diamesinae increases with the closer proximity to the terminus of melting glaciers [36]. Larocque et al. [35] found that "Pseudodiamesa and Diamesa were most abundant in alpine-tundra lakes above timberline, characterized by cold climatic conditions and low sedimentary organic content" in Swedish lakes. Protanypus has also been found in high elevation lakes in Canada and is associated with deep and cold lakes [35,80]. However, these taxa have not been identified in previous work done on modern chironomid distribution in the continental western USA [23,58,59,61,68]. Porinchu et al. [59] did find Pseudodiamesa in sediment collected from California from the interval corresponding to ages between 14,800 and 13,700 cal yr BP. No modern assemblages were comparable at that time, and the authors suggested that the presence of Pseudodiamesa indicated that the glacial meltwater was responsible for their deposition [59]. The presence and relative abundances of Diamesinae present in sediment collected from Rocky Mountain National Park suggest that this tribe may be used as a qualitative indicator of glacial meltwater and may assist historical reconstructions that use chironomids as a biological proxy for temperature.
DCA showed a strong relationship between glacial lakes and lakes only receiving melt from annual snowfall (Fig. 3). The first DCA axis possesses taxa associated with the colder lakes typical of glacial input in the negative range of DCA axis 1. Taxa affiliated with warmer temperatures are located to the right of the axis and are positive. DCA axis 2 represents the presence or absence of macrophytes in the system. Positive values are indicative of taxa typically affiliated with the presence of macrophytes (i.e., Psectrocladius, Paratanytarsus, Cladotanytarsus, and Tanytarsus). Negative values are affiliated with taxa that are typically found in the littoral zone of lakes or even small running streams such as Eukiefferiella, Diplocladius, Limnophyes, and Cricotopus/ Orthocladius. The top left quadrant of Fig. 3 contains some of the coldest stenotherms that have been noted in the literature. According to Brooks et al. [8], Abiskomyia only occurs in the coldest lakes of the arctic. Heterotrissocladius, Diamesa spp., and Pseudodiamesa are also noted as the coldest stenotherms present in assemblages [59,83]. Cony and Black Lakes are the only lakes that contain this assemblage. The taxa found in the bottom left quadrant are still indicative of cold water but are also affiliated with running water from streams or taxa more likely to be found in the littoral zone of lakes. Diplocladius, Limnophyes, Eukiefferiella, and Smittia are all uncommon in lake sediment and indicative of cold running water entering into the lake system. Surprisingly, two lakes that were thought to be only fed by year-of-snow, fall within this ordination space. Falcon Lake strongly falls within this zone. This lake was located in a rocky cirque and contained large snowfields that may mirror the action of glacially fed lakes. An unidentified rock glacier may also be feeding meltwater into this lake. Thunder Lake barely falls within this ordination space and may be indicative of a transitional lake. However, it may also suggest that this lake is receiving cold meltwater from somewhere higher in its catchment.

Relationship between environmental variables and chironomid communities
RDA was used to assess which environmental variable explains the most variance in species distribution [77] ( Table 3, Fig. 4). Only 2 of the 26 environmental variables were statistically significant with two others being close; SWT (p = 0.037), NO 3 + NO 2 -N (p = 0.049), boron (B) (p = 0.053), and %C (p = 0.057). The first axis of the RDA, consisting of SWT, boron, and %C, was identified as statistically significant (p = 0.022). Surface water temperature explained the majority of the variance present within the data (7.5%), followed by NO 3 + NO 2 -N (7.5%), B (7.4%), and %C (7.1%). The first RDA axis explained 9.8% in the variation of the distribution of chironomid taxa and was the only axis that was statistically significant (p = 0.019). Relationships between statistically significant variables were explored using Pearson's correlation coefficient. Only SWT and NO 3 + NO 2 -N were strongly and negatively correlated (r = −0.82, p = 0.007), and glacial meltwater is the driver that explains the relationship between the two variables. For this study, NO 3 + NO 2 -N was removed from the analysis as the correlation is too strong. However, this trend is well-documented in the western USA [9,50,73]. The story of atmospheric deposition of nitrogen in the western USA is a complicated one. While nitrogen has been collecting on glaciers for decades due to urbanization, fossil fuel consumption, and agriculture, biological activity is also active and contributes to this very complicated nitrogen story [16]. In addition, nitrogen-caused biological effects are related to contemporary atmospheric deposition to the east of the continental divide in Rocky Mountain National Park [87]. Elevated air temperatures in the latter part of the twentieth and into the twenty-first century have caused glaciers to recede, which has also introduced nitrogen into these systems [4]. Slemmons et al. [75] found that glacially fed lakes in the Rocky Mountains are 47 times higher in nitrogen than in snow-fed lakes. The Loch Vale Watershed is a long-term research site in Rocky Mountain National Park. Regional data collected from the Loch Vale Watershed site indicate that nitrate concentrations in alpine streams have increased by 50% since 2000 [4]. While SWT was the variable used in the analysis, this variable represents glacial meltwater contribution. The surface water temperatures of glacially fed lakes were ~2.62 °C colder than their paired lakes that only receive year-of-snow meltwater. The average difference in nitrogen was 66% higher in glacial lakes. While the relationship between atmospheric nitrogen deposition and algal communities is established [74,75], this is the first study to find a relationship between atmospherically deposited nitrogen and chironomids.
The relationship between chironomid communities and temperature is well established but poorly understood. Brundin (1949) noted cold stenotherms such as Heterotrissocladius spp. and Sergentia coracina present in late glacial sediment [60]. Quantitative paleotemperature reconstruction of SWTs was first performed by Walker et al. [81] on sediment collected from eastern Canada in the early 1990s. Air temperature models followed soon after. Lotter et al. [38] reconstructed air temperatures from a core collected from the Swiss Alps in the late 1990s [38]. However, the findings from this study indicate that SWT (i.e., glacial melt) is the environmental variable most responsible for the distribution of modern chironomid communities, whereas Mean July Air Temperature (MJAT) was not significant (p = 0.372). This finding illustrates that applying inference models to solely model air temperature may not always be appropriate when developing temperature reconstructions for fossil records as SWT and air temperature do not always covary. A need for a more in-depth understanding of how different temperatures affect the different lifecycles of chironomids is necessary to address future chironomid work. The results of this study suggest that active glacial activity present within a catchment will directly influence the chironomid communities present on the benthos. Future studies should acknowledge whether the lake under investigation is/or has ever been influenced by glacial meltwater. If the presence and/or absence is not known for the history of the lake, the presence of taxa from the tribe Diamesinae may act as qualitative indicator species downcore. Diamesinae may also indicate the presence of cold, flowing water into the system and would indicate that temperature reconstructions will produce colder conditions for SWT than were present for air temperature.
Initially, the relationship between boron and glacial activity was explored as a possible explanation for the distribution of boron in these lakes. However, the correlation between SWT and boron (r = −0.36, p-value = 0.345) as well as NO 3 + NO 2 -N and boron (r = 0.22, p-value = 0.566) was not statistically significant. Furthermore, two glacial lakes had boron levels below detection (Hutcheson and Odessa Lakes), and two snow-fed lakes have varying levels (Falcon Lake (0.05 ppm) and Thunder Lake (0.01 ppm)). The levels of boron are very low (0.00-0.05 ppm) and mirror minimum to median values for naturally occurring boron in surface water collected in British Columbia, Canada (0.01 ppm and 0.07 ppm, respectively) [48]. On average, levels of boron in US freshwater are about 0.10 ppm [10]. Maier and Knight [40] investigated the role of toxicity of waterborne sodium tetraborate on Chironomus decorus and found that growth rates were affected at 20 mg B/L, and acute toxicity occurred after 48 h at a level of 1376 mg B/L. However, the authors caution that aquatic macrophytes are much more susceptible to boron than macroinvertebrates and thus food dynamics for chironomids are more likely to be affected, which may explain the relationship between modern chironomid communities and the presence of boron in RMNP [40]. Other studies indicate higher uptake of boron in filamentous algae than chironomids [70]. Future research is needed to address these relationships in a more in-depth manner.
The presence of boron in lake water and its impacts on chironomid communities is poorly understood, and no studies currently exist that investigate this relationship to the author's knowledge. While boron may be a natural byproduct derived from weathering processes on sedimentary bedrock such as shales and coal deposits [3,48], the bedrock of all lakes sampled for this study is igneous diorites and granites interspersed with biotite gneisses and schists. However, boron may be used as an inorganic tracer of anthropogenic activity [3]. Boron is known as an indicator of wastewater and is a byproduct of nonchlorine bleach. The third pathway of boron deposition may result from the fly ash particles created from coal-fired power plants [3,48]. Due to the remoteness and elevation of the sample sites, this scenario is the most likely explanation for the presence of boron in high alpine lake water. The coal-fired Valmont Generating Station was located down the valley in the Boulder Creek Watershed from 1924 to 2017 and may be a possible source for the boron present in these systems. Trace amounts of boron could potentially be uplifted into the mountains by winds controlled by the summer monsoon. Additionally, the dominant westerly winds could be carrying boron derived from fertilizers to the west of the study site. Further research is needed to resolve the source of boron in this area. Carbon (%C) collected from bulk sediment is often used to understand the organic matter content within a lake [45]. Higher levels of organic carbon indicate that the lake is more productive, and larger sources of available food for chironomid larvae [19]. Many studies have found a strong statistical relationship between organic carbon and chironomids in Fennoscandia [52], Sweden [35], northwestern Canada [85], and New England of the USA [19]. Wilson and Gajewksi [85] argue that the large gradient of organic carbon collected for these studies captures a wide array (3-87%) that partially explain the distribution of chironomid communities in northern British Columbia and southwestern Yukon. The authors argue that other chironomid workers, such as Walker et al. [81], did not capture as full a gradient and were less likely to see this relationship. Our study only captures a gradient from 5.3 to 13.4%, and yet %C is an important environmental variable for understanding the modern distribution of chironomid communities in the alpine lakes of Rocky Mountain National Park. However, it should be noted that organic carbon content of bulk sediment is created from the interaction of primary productivity, allochthonous carbon, wind and wave action, sediment availability, distance from shore, light penetration, and nutrient availability [42,85]. For this reason, chironomid workers have been cautious in their interpretation of the relative importance of organic carbon as it is often highly correlated with depth and surface water temperature [38,52,85]. The relationship between surface water temperature and %C is correlated in this study (r = 0.66, p = 0.055). However, this relationship is not strong enough to warrant removing it from analysis as the statistical probability that the relationship between SWT and %C occurred by chance is more likely than the relationship evident between SWT and NO 3 + NO 2 -N.

Conclusion
The findings of this study indicate that glacial retreat is impacting the chironomid communities in the high elevation lakes located along the continental divide of the Colorado Rocky Mountains. Surface water temperature and NO 3 + NO 2 -N were extremely and strongly negatively correlated, indicating that glacial retreat is responsible for the greatest amount of explained variance (14.95%) from the model. Furthermore, limnological measures and the high presence of Chironomus and Sergentia suggest that early colonization of formerly oligotrophic to mesotrophic conditions is currently underway. However, cold stenotherms, such as Heterotrissocladius, are still present in high relative abundances suggesting these lakes are still affiliated with extremely cold conditions. The presence of taxa from the tribe Diamesinae (Diamesa, Pseudodiamesa, and Protanypus) are present in high numbers relative to the previous chironomid lacustrine studies and may indicate extremely cold and running water entering the lakes. These taxa may be useful as qualitative indicators of meltwater and may be useful for downcore paleotemperature chironomid-based reconstructions.
The findings from this study indicate that the high elevation lakes located in the remote lands of Rocky Mountain National Park are undergoing changes in trophic state, while a few are still maintaining the conditions that have been evident within these systems since the Pinedale glaciation. Almost all lakes in this study are no longer oligotrophic and are becoming more productive. The presence of boron in some lakes is also concerning as their presence indicates that anthropogenic activities are shaping these remote alpine ecosystems. This understanding will enable land managers for Rocky Mountain National Park to understand the current situation of water quality within the park.
The results from this study also inform our understanding of the processes that occur during the transition from glacial to interglacial stages in sediment. Many lakes that are studied for paleoclimatology are often found in remote locations and were formed by glacial activity. This study indicates that future work should endeavor to understand the glacial history within the lake catchment to refine midge-based temperature reconstructions. The presence of Diamesinae may suggest that warmer air temperatures were occurring as surface water temperatures were decreasing. Future research should explore the possibilities of combining reconstructions of surface water | https://doi.org/10.1007/s42452-021-04835-7 temperatures and air temperatures as drivers change within the system. The author suggests that the study sites presented here should be sampled after 2025 to record how these lakes have changed in the ten years since the study was undertaken. In addition, future paleoecological work should explore sediment that represents glacial retreat, such as those that followed the Pleistocene-Holocene transition.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.