Emergence of steeply stratified permafrost thaw ponds changes zooplankton ecology in subarctic freshwaters

ABSTRACT Climate change and associated permafrost thaw are creating new shallow waterbodies in vast regions of the circumpolar Arctic. These thaw ponds are characterized by high concentrations of colored dissolved organic matter originating from the degrading watershed, inducing a strong vertical thermal and oxygen (O2) stratification. We investigated the zooplankton community and biomass in eight subarctic thaw ponds and evaluated how zooplankton respond to this stratification. In a subset of three ponds, we further examined how other environmental variables, including essential fatty acids (EFA) concentration and phytoplankton, bacteria, and larval phantom midge Chaoborus biomass stratify and contribute to the vertical distribution of zooplankton in this increasingly common type of arctic freshwater system. The zooplankton community was extremely abundant in all ponds (up to 3,548 ind L−1) and dominated mainly by rotifers (35–93 percent of the biomass). Most zooplankton aggregated at the interface between the shallow well-oxygenated mixed surface layer and the deeper hypoxic but algal-rich stratified layer, and their distribution was affected by a combination of O2, Chaoborus, phytoplankton, and EFA that were supplied from opposite directions. Our findings show how water column stratification deeply affects the ecology of planktonic organisms in circumpolar freshwaters and indicate Arctic zooplankton species composition is expected to deeply change with the ongoing warming and browning.


Introduction
Frozen tundra soils represent a major component of circumpolar landscapes. They cover a considerable area of the Arctic (Zhang et al. 1999) and are one of the Earth's largest pools of organic carbon (Schuur et al. 2015). Ongoing climate warming in subarctic regions is causing this ice-rich permafrost to thaw rapidly and favor thermokarst processes that create shallow waterbodies, hereafter referred to as thaw ponds (Vonk et al. 2015). These freshwater systems have increased rapidly in number at high latitudes over the last decades, and these ponds now represent 25 percent of the area covered by lakes and ponds in the Arctic (Wauthy et al. 2018). They are strongly influenced by terrestrial inputs and act as recipients for large amounts of carbon and nutrients released from the degrading permafrost. As a result, water transparency and biogeochemical processes become highly altered (Wauthy et al. 2018), inducing a strong thermal stratification in thermokarst freshwaters and a shift from a polymictic regime to a dimictic regime (Breton, Vallières, and Laurion 2009). Consequently, as Arctic freshwaters become more dominated by thermally stratified waterbodies, there will be important implications for the Arctic freshwater species composition and ecology.
The growing influence of incoming terrigenous matter stimulates bacterial production ) and respiration (Roehm, Giesler, and Karlsson 2009) in these thaw ponds, making them highly net heterotrophic systems (Deshpande et al. 2015) and major greenhouse gas emitters to the atmosphere . Though the increasing amount of terrigenous nutrients also promotes pelagic primary producers (Przytulska et al. 2016;Wauthy and Rautio 2020), the high turbidity from these terrestrial inputs decreases light availability and induces a switch from benthic to pelagic primary production (Forsström, Roiha, and Rautio 2013). Furthermore, terrestrial inputs potentially reduce the nutritional quality of algae, decreasing in particular their essential fatty acid (EFA) content (Creed et al. 2018). In addition, the distribution of nutrients through the water column of thaw ponds is uneven; for example, Roiha, Laurion, and Rautio (2015) showed that total phosphorus levels were 6.5 times higher in bottom waters than at the surface. Chlorophyll displayed the same pattern, ranging from 9.5 µg L −1 at the surface to 133.4 µg L −1 in the bottom waters. However, rather than suggesting the presence of an abundant and active algal community in the deep hypolimnion of these deeper thermokarst freshwaters, this accumulation of chlorophyll has been interpreted as a consequence of slow decay of dead algae in the absence of light, low temperatures, and anoxic conditions as well as the presence of photosynthetic green sulfur bacteria that synthesize chlorophyll b and bacteriochlorophyll (Rossi, Laurion, and Lovejoy 2013;Roiha, Laurion, and Rautio 2015). Prokaryotes and heterotrophic nanoflagellates may also have a strong vertical gradient (Roiha, Laurion, and Rautio 2015;Deshpande et al. 2016), and the substantial bacterial production coupled to the steep thermal stratification can induce severe hypoxia in the bottom waters of thaw ponds (Deshpande et al. 2015). The result has a marked effect on pluricellular oxygen-dependent organisms living in thermokarst systems.
Despite these major changes in circumpolar freshwaters due to permafrost thaw, the repercussions on the higher trophic levels have just been addressed superficially. Only two studies have investigated the pelagic consumers in thaw ponds, emphasizing the high abundance of zooplankton and the dominance of rotifers (Nigamatzyanova and Frolova 2016;Bégin and Vincent 2017). In this study, we provide further details about zooplankton in thaw ponds by examining in detail the vertical distribution of zooplankton in thaw ponds. The main objectives were to (1) investigate whether and how the zooplankton community responds to the thermal and oxygen (O 2 ) stratification in subarctic thaw ponds, which is an unconventional physical characteristic for subarctic shallow ponds; (2) examine the effect of stratification on other environmental variables, including phytoplankton biomass and EFA concentration, as well as bacteria and larval phantom midge Chaoborus biomass; and (3) determine the main variables controlling zooplankton vertical distribution through the water column. Given the anticipated strong thermal structure and hypoxia observed in the hypolimnion of most sampled thaw ponds (Deshpande et al. 2015), we hypothesized a pronounced stratification of zooplankton. This stratification was expected to be the result of a compromise between a good quality algal diet, O 2 accessibility, and top-down predation.

Study sites
Fieldwork was conducted from 4 to 8 August 2013 in two permafrost thaw-affected valleys located in an area of sporadic permafrost near the village of Kuujjuarapik-Whapmagoostui, northern Quebec, Canada (55°17′ N, 77°47′ W). The Kwakwatanikapistikw River valley (KWK) is filled with numerous thaw ponds formed in clay-silt bed depressions created by the thawing of mineral palsas (Bouchard et al. 2011). Their watershed is mainly colonized by dense shrubs (Salix spp., Alnus crispa, Myrica gale) and sparse trees (e.g., Picea mariana, P. glauca, Betula glandulosa), with some patches of mosses (Sphagnum spp.) and macrophytes (Carex spp.; Bouchard et al. 2011). The Sasapimakwananisikw River valley (SAS) is a peatland that is covered predominantly by mosses (Sphagnum spp.) and macrophytes (Carex spp.) and surrounded by thermokarst waterbodies and collapsing organic palsas (Arlen-Pouliot and Bhiry 2005). Pictures of the study sites are available in the Supporting Information Figure S1. Further details about KWK and SAS are provided in Bouchard et al. (2014) and Arlen-Pouliot and Bhiry (2005). Most of the thaw ponds in KWK and SAS are small (mean surface area 183 ± 46 m 2 ) and shallow (mean depth 2.3 ± 0.3 m) freshwaters displaying high nutrient and suspended organic and inorganic matter concentrations, especially in the bottom waters (Table 1).

Sample collection
In total, eight ponds were sampled, all fishless (KWK1, KWK2, KWK6, KWK12, SAS2A, SAS2B, SAS2C, SAS2D). We measured profiles of temperature, dissolved O 2 concentrations, and photosynthetically active radiation through the water column of each pond using a YSI 550A probe (YSI Inc., Yellow Springs, Ohio) and a PUV-2500 profiler radiometer (Biospherical Instruments, San Diego, California). Based on the thermal stratification and aiming for samples from the mixed and stratified layers, we collected water from four to five depths in 40-cm to 1-m intervals between 0 m (just below the surface) and about 0.2 m above the bottom (Figure 1), using a 2-L transparent Plexiglas water sampler (diameter of 10 cm, length of 30 cm, Limnos Ldt., Turku, Finland).
To be able to sample in fine vertical resolution, the sampler was used in a horizontal position as in Rautio, Korhola, and Zellmer (2003) when collecting water from the surface depths (0, 0.4, and 0.8 m depth). All ponds were sampled for chlorophyll a and a subset of three ponds for fatty acids, phytoplankton, and bacteria. Zooplankton and Chaoborus phantom midge larvae were sampled with same sampling bottle from a 6-L water sample (triplicates per depth) that was concentrated using a 50-µm plankton net. The organisms were preserved in formaldehyde (4 percent final concentration) until needed for counting. The sampling took place between 12:00 and 16:00 to standardize the effect of potential diel vertical migration on plankton depth distribution.
In the lab, 200-mL duplicate subsamples of water from each pond were filtered through GF/F glass fiber filters (0.7 μm, 47 mm, Whatman plc, Maidstone, United Kingdom) and stored at −80°C until extraction in 95 percent ethanol to determine chlorophyll a concentration following the protocol of Nusch (1980). For three ponds (KWK6, KWK12, and SAS2A), water was also filtered to investigate bacterial and phytoplankton biomass and fatty acid concentrations. Prior to filtration, bacterioplankton were preserved by adding prefiltered (0.2 µm) glutaraldehyde (2 percent final concentration) to a 10-mL sample (triplicates) and stored at 4°C for 12 hours. Each sample was then stained with 4′,6-diamidino-2-phenylindole (5 µg mL −1 final concentration) for 5 minutes and filtered through a nucleopore black polycarbonate membrane (0.22 μm, 25 mm, Whatman plc) placed on the top of a GF/F glass fiber filter (0.7 μm, 25 mm, Whatman plc). The nucleopore membrane was then mounted on a microscope slide covered with nonfluorescent immersion oil and a coverslip and kept at −20°C until counting. Phytoplankton were preserved with acid Lugol's solution (Sournia and Caspers 1980) and stored at 4°C until counting. Water samples for fatty acids (150-500 mL) were filtered onto precombusted and preweighted GF/F glass fiber filters (0.7 μm, 47 mm, Whatman plc). We visually inspected the filters and removed zooplankton and larger detritus before storing the filters at −80°C.

Zooplankton community and biomass
Zooplankton were identified and enumerated using the Utermöhl sedimentation method (Utermöhl 1958). Stemberger (1979), Hebert (1995), and Haney et al. (2013) were used as taxonomic resources to identify specimens to morphological species. Counts were made under an Axio Observer.A1 inverted microscope (Carl Zeiss, Oberkochen, Germany) at 100× to 400× magnification. At least 400 specimens were counted for each sample. To estimate the mean zooplankton biomass, we measured the body width and length of ten individuals per morphological species using microphotographs captured by an Axiocam ERc 5S coupled to the AxioVision Rel. 4.8.2. software (Carl Zeiss). Rotifer biovolumes were estimated as per McCauley (1984). A 1:1 ratio was considered between the biovolume (mm 3 ) and the wet weight (mg). This was converted to dry weight using a ratio of 0.1 as in Pace and Orcutt (1981). We calculated the biomass of the cladoceran Bosmina spp. and cyclopoid copepod nauplii according to Rosen (1981), and we used the regression equations of Malley et al. (1989) to estimate the biomass of Holopedium gibberum and cyclopoid copepodites CI-CV. The biomass of the cladoceran Daphnia pulex s.l. was estimated as in Dumont, Van de Velde, and Dumont (1975).

Bacteria, Chaoborus and phytoplankton biomass
We analyzed bacterioplankton samples by epifluorescence microscopy using an Axio Observer.A1 inverted microscope (Carl Zeiss) at 1,000× magnification under UV light. The bacteria biomass was calculated by multiplying the bacterial abundance with the average bacteria cellular carbon content (18.1 fg C cell −1 ) as in Roiha, Laurion, and Rautio (2015). We counted the dipteran larvae Chaoborus sp. from zooplankton samples and measured their body length under a SteREO Discovery.V12 microscope (Carl Zeiss). We then Table 1. Limnological properties (mean ± standard error) of subarctic thaw ponds sampled in this study in Kwakwatanikapistikw (KWK) and Sasapimakwananisikw (SAS) Rivers valleys. Note. The table displays variables sampled from the surface and the bottom of the water column, including dissolved organic carbon (DOC), particulate organic carbon (POC), total suspended solids (TSS), total nitrogen (TN), total phosphorus (TP), soluble reactive phosphorus (SRP), pH, and conductivity (Cond.). Source. Data were retrieved from Breton, Vallières, and Laurion (2009)  calculated their biomass following the regression equation of Mason (1977). Phytoplankton were identified and counted following Utermöhl (1958) and calculated biomass from biovolumes as per Forsström et al. (2015).

Fatty acid analyses
Seston fatty acids were extracted from freeze-dried samples as per Mariash et al. (2011) using a redesigned extraction method from Bligh and Dyer (1959) and were methylated following an adapted protocol from Lepage and Roy (1984). First, a mixture of methanol-toluene and acetyl chloride (4:1:0.125) was added to the samples with internal standard (nonadecanoic acid C19:0, N5252, available from Sigma-Aldrich, Saint Louis, Missouri). After centrifugation, the samples were incubated at 90°C for 20 minutes. Transesterified fatty acids were extracted with hexane and were scanned through a 7890A gas chromatographer (Agilent, Santa Clara, California) coupled to a 5975C mass spectrometer with a triple-axis detector (Agilent) and a J&W DB-23 column (Agilent). Identification and quantification were done using the calibration curves that were based on different concentrations of a standard mix (thirty-seven fatty acids components; Supelco 47885-0, Sigma-Aldrich). The eicosapentaenoic acid (EPA) and docosahexaenoic acid (DHA) are polyunsaturated omega-3 fatty acids synthesized almost exclusively by algal primary producers, and they are required for the normal growth and reproduction of most heterotrophic organisms (Taipale et al. 2014). Because zooplankton are unable to produce EPA and DHA de novo, they need to obtain these fatty acids from their diet (Parrish 2009). Consequently, in this study, the concentration of EPA and DHA were summed and designated as EFA.

Data analyses
We statistically estimated the depth of greatest change in temperature and dissolved O 2 for each pond using these two variables as a bivariate matrix. Data were normalized and cluster analyses (with group average cluster mode), followed by a similarity profile analysis (SIMPROF) routine, were carried out to identify the bivariate structure in the data against the depth gradient. A significant change in structure was considered as a boundary layer in temperature and dissolved O 2 (hereafter referred to as cline), separating a mixed layer in the surface and a stratified layer in deeper waters. The weighted mean depth relative to this cline was calculated for each zooplankton group (Rotifera, Copepoda, Cladocera) in each sampled pond using the following adapted equation of Worthington (1931): where D is the depth of the cline (m) in the sampled pond, d i is the sample depth (m), and B i is the biomass of each zooplankton group (µg L −1 ) measured in d i . Because the sampling depth intervals were not equal over the water column, l i was introduced to the equation and calculated as the distance (m) from half the distance to the previous sampling depth (or surface) to half the distance to the next sampling depth (or bottom). Taxon differences in weighted mean depths were analyzed by a Kruskal-Wallis rank test, with zooplankton taxa as independent variables, using R v3.5.1. (R Development Core Team 2018).
We investigated the relationships between the vertical distribution of zooplankton-total zooplankton, rotifer, and copepod biomass-and the potential explanatory variables of EFA concentration, chlorophyll a concentration and phytoplankton biomass (algal variables), Chaoborus and bacteria biomass (nonalgal variables), and temperature, dissolved O 2 concentrations, and photosynthetically active radiation (physicochemical variable) in three ponds (KWK6, KWK12, SAS2A) using distance-based linear models (DistLM). Due to their absence from KWK6 and KWK12, cladocerans were not tested separately in these analyses. Zooplankton biomass was square root transformed, and resemblance matrices were generated using Bray-Curtis similarity. Explanatory variables were square root transformed and normalized. The collinearity between explanatory variables was also tested via a draftsman plot, and variables that were highly correlated (r > |0.80|) were removed from the analysis. Thus, chlorophyll a concentration was omitted from the model due to its high correlation with phytoplankton biomass (r = 0.86), and temperature and photosynthetically active radiation were removed due to their high correlation with dissolved O 2 concentrations (r = 0.83 and 0.81, respectively). Due to the small numbers of samples (thirteen), we used the corrected Akaike information criterion (AIC c ) to rank models for all possible combinations of explanatory variables, with lower AIC c values and higher AIC c weights indicating better model performance (Burnham and Anderson 2002). DistLM analyses were run using PRIMER v6.1.16. with the add-on package PERMANOVA+ v1.0.6. (PRIMER-E Ltd., Plymouth, United Kingdom). The number of permutations was fixed at 9,999.

Water column stratification
Despite their shallow depths (ranging from 1.9 to 3.2 m, median depth of 2.7 m), all sampled thaw ponds showed strong physicochemical stratification. Temperatures decreased 7.9 ± 0.7°C, and dissolved O 2 concentrations dropped 7.7 ± 0.5 mg L −1 between the surface and bottom waters (Figure 1). Photosynthetically active radiation decreased rapidly with depth in all ponds, especially in the SAS ponds, where the 1 percent surface radiation depth averaged 1.06 ± 0.05 m (Figure 1). In all ponds, O 2 was under saturated even at the surface. In the clearer ponds (KWK1, KWK2, KWK6, KWK12), the first half of the water column was thermally homogenous and well oxygenated (9.1 ± 0.2 mg L −1 ), whereas in the darker ponds (SAS2A, SAS2B, SAS2C, SAS2D) the mixed layer was shallower and displayed O 2 concentrations ranging between 2.3 and 7.6 mg L −1 (Figure 1). The depth of the cline separating the mixing and stratified layers varied between ponds, ranging from 0.3 to 2.3 m. All ponds were hypoxic at the bottom (0.4 ± 0.1 mg L −1 ).
Chlorophyll a concentrations varied widely among the thaw ponds (Figure 1), both at the surface (2.3 to 19.9 µg L −1 ) and in deeper water (16.2 to 196.5 µg L −1 ), but had the same profile pattern. In all ponds, chlorophyll a concentrations increased with depth and were from 2.2 to 45.1 times higher at the bottom than at the surface.

Zooplankton community and biomass
We identified thirty-eight zooplankton taxa from all ponds combined, with fourteen to twenty-five taxa observed per pond (Supporting Information Table  S1). Rotifers were the most diverse group, with thirtyfour taxa. We observed only three Cladocera and one Copepoda taxa in all ponds. In ponds KWK6 and KWK12, we observed no cladocerans. The total abundance of zooplankton averaged 442 to 3,548 ind L −1 , and rotifers represented 97.6 ± 0.6 percent of this total abundance (Supporting Information Table S1).
The average zooplankton biomass per pond was 42 to 529 µg L −1 (Table 2). Rotifers dominated the zooplankton assemblages, representing 34.6 to 92.6 percent of the total biomass. The exception was pond SAS2B, which was dominated by D. pulex s.l. Eight taxa represented more than 95 percent of the rotifer biomass: Polyarthra spp.

Vertical distribution of zooplankton and environmental variables
The weighted mean depth of each zooplankton group indicated a similar pattern of vertical alignment within the water column (Kruskal-Wallis, p = .19), displaying a preference for the interface between the mixed and stratified layers (Figure 2). Overall, in most ponds, zooplankton were located within 0.5 m of the cline that determined the greatest change in temperature and O 2 .
In the three ponds (KWK6, KWK12, and SAS2A) where we examined zooplankton distribution in greater detail and associated the zooplankton with the physicochemical and biological variables, we noted profiles that followed the same overall pattern (Figure 3). The highest biomass of zooplankton was found in the mid-water column close to the cline. The very surface displayed a low biomass of zooplankton, and the biomass also declined when O 2 concentration was very low. The maximum biomass of zooplankton, however, varied among the ponds and was nearly five times higher in SAS2A than in KWK6.
In these three ponds, phytoplankton biomass varied from 0.1 to 70.2 mg L −1 but displayed similar stratification patterns, with the lowest values in the mixed layers and the highest values in the stratified layers ( Figure 3). EFA had very low concentrations at the surface, then increased rapidly to reach maximum values at the cline, and finally decreased in the stratified layers. Bacteria and Chaoborus biomasses did not present any distinct generalized vertical distribution patterns among the ponds. No Chaoborus were found in SAS2A.
The vertical distribution of zooplankton had a clear, positive relationship with O 2 concentrations, as was reflected in the model fit values of R 2 = 0.29 for total zooplankton, R 2 = 0.25 for rotifers, and R 2 = 0.16 for copepods (Table 3). Iterating Chaoborus, phytoplankton, and EFA into the models improved the fit, at R 2 = 0.36 for total zooplankton, 0.28 for rotifers, and 0.29 for copepods. For total zooplankton, the best model was based on the relative likelihood of O 2 alone (depicted by the AIC c weight of 0.25), whereas the second (Chaoborus) and third (O 2 + phytoplankton) best models had respective AIC c weights of 0.11 and 0.09. The vertical distribution of rotifers was best explained by the models that included O 2 , phytoplankton, and O 2 + phytoplankton; however in these models, the O 2 -alone model explained the zooplankton distribution 3.4 times better than the next model based on phytoplankton biomass (AIC c weights of 0.31 and 0.09, respectively) as well as explaining most of the R 2 . The other variables had a much smaller effect. For the copepods, O 2 and Chaoborus biomass alone constituted the two best models, having similar AIC c and AIC c weights; however, EFA was an important variable by doubling the R 2 when integrated (Table 3). When considered alone, only O 2 was significant (DistLM marginal test, p < .05), explaining 29 percent of the overall variation of total zooplankton (Supporting Information Table S2).

Discussion
We found large variation in temperature and O 2 of thaw ponds between the mixed surface layer and the deeper strata. These were consistent across the eight ponds studied. Such physicochemical stratification is caused by permafrost thaw and comes with important implications for the organisms living in these systems. We demonstrated how in all ponds this stratification explained the accumulation of zooplankton to the interface between the warmer, oxygen-rich layer and the hypoxic but resource-rich layer. We will further discuss how the overall zooplankton community composition was determined by the O 2 concentration that likely favors rotifers over crustacean zooplankton.

Water column stratification
Consistent with previous studies (e.g., Breton, Vallières, and Laurion 2009;Deshpande et al. 2015), we observed a steep thermal stratification within the pelagic zone of all sampled thaw ponds (Figure 1). Such a strong temperature gradient in shallow freshwaters can be explained by the high concentration of chromophoric dissolved organic matter (Roiha, Laurion, and Rautio 2015;Wauthy et al. 2018) and total suspended solids (Watanabe et al. 2011;Crevecoeur et al. 2015) that considerably decrease light penetration through the water column due to diffusion and absorption (Figure 1). Given this strong thermal stratification, gas exchange is limited between the stratified layer and atmosphere during the summer (Deshpande et al. 2015), thereby contributing to the O 2 depletion in bottom waters (Figure 1). In all ponds, the bottom waters were hypoxic, and in some ponds the hypoxic layer extended upward, leaving only the very surface oxygenated.
In addition to marking the distinct vertical change in temperature and O 2 , the chlorophyll a profile was also very different between the top and bottom layers of the water column (Figure 1). The already high concentration of chlorophyll a in the surface waters in comparison to clearwater Arctic ponds  increased substantially toward the bottom. As in earlier studies (Roiha, Laurion, and Rautio 2015), we suggest that this high concentration of chlorophyll a is a combination of cold and hypoxic waters that slow down the degradation of sedimenting algae and the presence of photosynthetic green sulfur bacteria that synthesize chlorophyll b and bacteriochlorophyll (Rossi, Laurion, and Lovejoy 2013;Przytulska et al. 2016). These compounds can interfere with algal chlorophyll a measurements in fluorometric methods (Tolstoy and Toth 1980;Coveney 1982). The fact that the chlorophyll a concentrations were lowest in KWK6, the only pond where most of the water column was well oxygenated, supports this explanation. Furthermore, the high dissolved organic matter concentration contributed to light limitation (Figure 1), which presumably inhibited primary production in deeper depths. However, even inactive, the high chlorophyll a values indicate a high availability of putative resource biomass for zooplankton in bottom waters.

Zooplankton communities in thaw ponds
Possibly due to the limited habitat size of welloxygenated water, the thaw pond zooplankton Table 2. Mean water column biomass (µg L −1 ) of zooplankton, including Rotifera, Copepoda, and Cladocera in the sampled subarctic thaw ponds. community composition favored a high abundance of rotifers. These organisms are known to be tolerant of low O 2 concentrations (e.g., Bērziņš and Pejler 1989;Mikschi 1989), although they, too, are absent in thaw ponds under winter ice (Rautio, unpublished) when the entire water column is anoxic (Matveev, Laurion, and Vincent 2019). The dominance of rotifers compared to planktonic crustaceans in our study is consistent with previous studies from thermokarst sites (Nigamatzyanova and Frolova 2016;Bégin and Vincent 2017). Moreover, despite their small size, rotifers contributed substantially to the overall zooplankton biomass (Table 2). This supports the observations of Bégin and Vincent (2017) highlighting the rotifer dominance in thermokarst systems and is in contrast with arctic (Rublee 1992), subarctic (Pinel-Alloul et al. 1982), and temperate lakes (Herzig 1987) with higher O 2 concentrations.
The dominance of rotifers can further be explained by their pelagic feeding mode and their ability to graze very efficaciously on diverse sources of diet such as algae, bacteria, and even particulate organic matter (Arndt 1993). These resources are abundant in thaw ponds (Roiha, Laurion, and Rautio 2015) that, unlike clearwater arctic ponds and lakes ), can be considered mesotrophic based on their nutrient levels (Table 1) and high chlorophyll a concentration (Figure 1). Previous studies have highlighted the ability of rotifers to gain an initial advantage when nutrients and phytoplankton biomass increase (Marty, Pinel-Alloul, and Carrias 2002;Jalal, Pinel-Alloul, and Méthot 2005). It is, however, important to note that though new thaw ponds are constantly been formed in permafrost thaw regions, the studied ponds are estimated to be older and were formed more than fifty years ago (Fillion, Bhiry, and Touazi 2014). The zooplankton communities inhabiting these waters should therefore be considered established. Rotifers could have formed the initial community but continued to dominate as the ponds aged, most likely due to the constant supply of nutrients and dissolved organic matter from the thawing catchment, which provides a physicochemical setting that is favorable to these zooplankton. In addition, the short life span of rotifers and high reproductive rates make these r-strategy zooplankton particularly well adapted to perturbed environments (Jalal, Pinel-Alloul, and Méthot 2005). Although the suppression of rotifers by large Daphnia populations through exploitative competition for food resources and mechanical interference is well documented (Gilbert 1988), rotifers displayed a high biomass in our study sites, despite the presence of  Daphnia. This suggests that there was no significant influence of cladocerans on rotifer biomass, possibly because the Daphnia were first controlled by the overall low O 2 concentration in these ponds.
On average, we identified seventeen rotifer taxa per thaw pond (Supporting Information Table S1). This rotifer diversity is ten taxa higher than that found in a previous study from thaw ponds in the same region (Bégin and Vincent 2017). The eight rotifer taxa that dominated the thermokarst zooplankton communities are common throughout northern Quebec (Pinel-Alloul et al. 1982;Bégin and Vincent 2017). These taxa can rely on various food sources, including phytoplankton and other species of rotifers but also bacteria and heterotrophic flagellates (Pourriot 1977;Arndt 1993). However, the taxon Polyarthra spp., which represented almost 73 percent of the rotifer biomass, are selective feeders of larger particles, such as pelagic algae (Pourriot 1977).
In contrast with the rotifers, zooplankton crustaceans displayed a very low diversity in the sampled thaw ponds. We found only three members of Cladocera, a diversity much less than that usually observed in clearwater ponds and lakes of subarctic Canada (Pinel-Alloul et al. 1982;Paterson et al. 2014). The most dominant cladoceran we observed was D. pulex s.l., a cosmopolitan species able to feed on algal (Pennington 1941) and bacterial (Peterson, Hobbie, and Haney 1978) food sources. For the subclass of Copepoda, our taxonomic resolution was at the order level due to the absence of adult copepods in the samples, limiting a more precise identification. However, all encountered copepod larvae (nauplii) and juveniles (copepodites) were cyclopoids. Microcyclops was the only genus of the Cyclopoida order observed in KWK and SAS sites in a previous study (Bégin and Vincent 2017), whereas eleven species of copepods are known to exist in the clear ponds and lakes of the region (Swadling et al. 2001).
Because all of our samples were collected from the center of the ponds and during the day, we cannot exclude zooplankton horizontal movement as a contribution to the observed species composition. In particular, cladocerans can migrate horizontally between offshore and littoral zones, especially in small and shallow freshwater systems (Burks et al. 2002). Macrophytes along the shoreline of all of our thaw ponds serve as potential refuges for zooplankton from Chaoborus predation pressure that can limit cladoceran abundance in the center of the ponds during the day (Kuczyńska-Kippen 2001). This scenario of horizontal migration, however, is unlikely to explain the complete absence of common cladoceran taxa, such as Bosmina, in some ponds and the generally low abundance of copepods and cladocerans as well as the overwhelmingly dominant abundance of rotifers in these small fishless ponds.

Vertical distribution of zooplankton
The weighted mean depth results indicated that all zooplankton taxa stayed preferentially at the interface between the mixed surface layer and hypoxic deeper layer (Figure 2). If O 2 was the only environmental variable controlling the vertical distribution of zooplankton, we would have expected zooplankton to remain in the warmer and best-oxygenated waters at the surface of the mixed layer. However, zooplankton appeared drawn to the depths where phytoplankton, and presumably their diet, was more abundant.
The more detailed analyses in KWK6, KWK12, and SAS2A highlighted the environmental variables potentially driving vertical distribution of zooplankton. Based on the DistLM models, dissolved O 2 concentration was the most important environmental predictor of total zooplankton, rotifer, and copepod distribution (Table  3). Given the presence of cladocerans solely at the oxygenated surface (Figure 3), dissolved O 2 would seem to be the main driver of cladoceran vertical distribution as well. Temperature was not included in the linear models due to its high correlation with O 2 , which prevented us from testing its unique contribution to zooplankton distribution. The importance of this variable has been tested for zooplankton growth in earlier thaw pond studies (Przytulska et al. 2015), and it has been shown that the growth response was complex and depended on species-specific adaptations and interaction with food availability, but in general the thaw pond zooplankton grew best at 18°C. We therefore assume that the warmer mixed layer (15°C) would stimulate growth more than the cold bottom waters (<10°C) and that temperature contributed in making the mixed layer favorable to zooplankton. Our sampling was restricted to the daytime, and zooplankton are known to migrate diurnally throughout the water column (Hutchinson 1967). The patterns of diel vertical distribution displayed by rotifers and copepods are, however, known to be rather consistent in these ponds (Bégin and Vincent 2017), which is expected in a situation where O 2 is the driving environmental variable; its profile does not change in diel bases but stays constant (Deshpande et al. 2015).
The models also indicated that top-down predation applied by Chaoborus explained an appreciable proportion of variation. Thus, O 2 concentration and predator avoidance drove much of the vertical distribution of the total zooplankton community. In the absence of fish and adult copepods, Chaoborus larvae are expected to be the main predators of zooplankton, for both microzooplankton and cladocerans (Lewis 1977;MacKay et al. 1990). Chaoborus biomass varied greatly among ponds, ranging from 0 to 1,360 µg L −1 . However, our sampling method was not ideal for capturing Chaoborus, and their abundance may have been underestimated due to too low sample volume. It could also be that the low number of Chaoborus in our samples was due to their vertical and horizontal migration, which is well documented for phantom midge larvae (Voss and Mumm 1999). Bégin and Vincent (2017) emphasized earlier the low presence of Chaoborus in SAS2A samples collected at noon, whereas they observed higher numbers in samples collected around midnight. This pattern reflected the trait of Chaoborus to hide in the sediments and littoral zones during the day and migrate back to the pelagic zone during the night (Hare and Carter 1986).
The models also included the algal food sources (phytoplankton and EFA and indirectly chlorophyll a, which was correlated with phytoplankton biomass) among the best models, reflecting the influence of diet on the vertical distribution of zooplankton. Phytoplankton biomass was highest in the bottom waters of the thaw ponds (Figure 3), most probably due to the same reasons as the high chlorophyll concentration; that is, the slow decay rates of sedimenting phytoplankton cells. The degrading nature of phytoplankton in deeper waters is also supported by the concentration of EFA in seston. EFA had a low concentration in the surface and increased only to a certain depth that was near the interface between the illuminated surface and nutrient-rich deeper layers, before decreasing close to the dark bottom ( Figure 3). The best models did not include bacterial biomass as an explanatory variable for any of the taxonomic groups. Bacteria showed high abundance in all ponds, displaying biomass typical of thaw ponds (Roiha, Laurion, and Rautio 2015), and they were probably not limiting zooplankton diet in any depth. Further, although rotifers and cladocerans can feed on bacteria, the higher quality algal diet is likely more important in determining the vertical position of zooplankton in the water column, as was suggested by the statistically significant contribution of EFA to copepod distribution.

Conclusions
Climate change and the associated thawing permafrost are modifying northern landscapes and drastically altering hydrological regimes, biogeochemical processes, and biological production in arctic and subarctic lakes and ponds (Wrona et al. 2016;Wauthy et al. 2018). Presently, rotifers appear to benefit from these rapid changes, with subarctic thaw ponds containing a high biomass of these pelagic micrograzers. The strong thermal structure of thermokarst lakes and ponds, which is unusual compared to the normally polymictic subarctic freshwater, leads to a pronounced stratification of O 2 that drives the vertical distribution of zooplankton and possibly the overall species composition as well by selecting and organizing species according to their tolerance to hypoxic conditions. Furthermore, permafrost thaw releases organic matter and nutrients from the eroding watershed, stimulating phytoplankton and pelagic primary production despite the light limitations (Vonk et al. 2015). Grazer species that are able to cope with hypoxia can benefit from the high phytoplankton production in these nutrient-rich systems. Future warming (Wrona et al. 2016) and browning (Wauthy et al. 2018) of circumpolar freshwaters are expected to enhance the thermal stability of thaw ponds and therefore increase the potential of anoxia throughout the water column (Deshpande et al. 2015). This scenario adds much uncertainty in terms of the future response of arctic and subarctic zooplankton communities to climate warming and degrading permafrost.

Disclosure statement
No potential conflict of interest was reported by the authors.