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 (O₂) 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⁻¹) 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 O₂, 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.

KEYWORDS:

• Thaw ponds
• zooplankton
• stratification
• hypoxia
• top-down predation
• phytoplankton

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 (Deshpande et al. 2016) 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 (Matveev et al. 2016). 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⁻¹ at the surface to 133.4 µg L⁻¹ 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₂) 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₂ accessibility, and top-down predation.

Methods

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²) 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).

Table 1. Limnological properties (mean ± standard error) of subarctic thaw ponds sampled in this study in Kwakwatanikapistikw (KWK) and Sasapimakwananisikw (SAS) Rivers valleys

Strata	DOC (mg L^−1)	POC (mg L^−1)	TSS (mg L^−1)	TN (mg L^−1)	TP (µg L^−1)	SRP (µg L^−1)	pH	Cond. (µS cm^−1)
Surface	12.6 ± 1.2	1.9 ± 0.4	12.2 ± 1.8	0.4 ± 0.0	46.7 ± 4.0	2.5 ± 0.4	6.3 ± 0.2	59.4 ± 7.7
Bottom	11.2 ± 1.3	5.8 ± 1.2	42.4 ± 10.5	1.0 ± 0.2	228.3 ± 51.5	17.2 ± 8.2	6.1 ± 0.1	200.2 ± 21.6

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), Laurion et al. (2010), Bouchard et al. (2012), Crevecoeur et al. (2015), Deshpande et al. (2015, 2016), Roiha, Laurion, and Rautio (2015), and Wauthy et al. (2018).

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₂ 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.

Figure 1. Vertical profiles of temperature (°C), dissolved oxygen (O₂, mg L⁻¹), photosynthetically active radiation (PAR) irradiance (%), and chlorophyll a (µg L⁻¹) in the sampled subarctic thaw ponds. The horizontal dashed line defines the cline based on temperature and dissolved oxygen, between the mixed and stratified layers

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⁻¹ 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³) 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⁻¹) 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 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₂ 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₂ (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):$Weightedmeandepth=D-\text{ left({{{}\sum \text{^ }^}}\left(d_i{B}_i{l}_i\right)/\text{{{}\sum \text{^ }^}}(B_i{l}_i)\text{} right)}$,

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⁻¹) 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₂ 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₂ 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.

Results

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₂ concentrations dropped 7.7 ± 0.5 mg L⁻¹ 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₂ 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⁻¹), whereas in the darker ponds (SAS2A, SAS2B, SAS2C, SAS2D) the mixed layer was shallower and displayed O₂ concentrations ranging between 2.3 and 7.6 mg L⁻¹ (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⁻¹).

Chlorophyll a concentrations varied widely among the thaw ponds (Figure 1), both at the surface (2.3 to 19.9 µg L⁻¹) and in deeper water (16.2 to 196.5 µg L⁻¹), 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 thirty-four 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⁻¹, 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⁻¹ (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. (73.0 percent), Keratella testudo (5.3 percent), Trichocerca spp. (5.0 percent), Ascomorpha spp. (4.0 percent), K. cochlearis (3.4 percent), K. tecta (3.3 percent), Asplanchna spp. (2.6 percent), and Kellicottia longispina (1.7 percent). In the subclass of Copepoda, we only observed Cyclopoida, and most specimens were at the nauplii stage. D. pulex s.l. was the most dominant cladoceran taxon in terms of biomass, except in pond KWK2, where H. gibberum was the main species.

Table 2. Mean water column biomass (µg L⁻¹) of zooplankton, including Rotifera, Copepoda, and Cladocera in the sampled subarctic thaw ponds

Taxa	KWK1	KWK2	KWK6	KWK12	SAS2A	SAS2B	SAS2C	SAS2D
Rotifera	58.34	56.34	38.74	37.92	127.77	19.89	271.32	87.82
Ascomorpha spp.	0.05	4.09	0.20	4.75	0.54	0.12	27.49	0.39
Asplanchna spp.	2.60	3.93	0.70	2.18	2.02	0.00	0.00	0.00
Kellicottia longispina	1.56	1.85	0.35	1.96	1.79	0.00	0.00	0.43
Keratella cochlearis	0.73	1.00	2.96	5.76	0.63	0.06	0.15	0.18
Keratella tecta	2.61	7.31	0.20	2.38	0.00	0.16	2.54	0.10
Keratella testudo	1.15	0.99	0.02	0.00	8.35	1.25	6.61	20.34
Polyarthra spp.	49.05	33.68	33.86	6.89	111.87	17.61	229.76	64.78
Trichocerca spp.	0.19	2.03	0.00	12.42	0.91	0.18	2.95	0.42
Copepoda	17.50	22.09	3.10	28.29	46.43	14.42	195.95	18.96
Cyclopoida nauplii	15.42	20.52	3.10	28.29	44.87	14.42	189.27	16.80
Cyclopoida copepodites	2.08	1.57	0.00	0.00	1.56	0.00	6.68	2.16
Cladocera	0.73	7.80	0.00	0.00	14.18	23.18	62.01	16.71
Bosmina sp.	0.00	0.00	0.00	0.00	0.00	0.00	23.48	3.06
Daphnia pulex	0.73	1.48	0.00	0.00	12.63	23.18	38.54	11.47
Holopedium gibberum	0.00	6.32	0.00	0.00	1.55	0.00	0.00	2.18
Total zooplankton	76.57	86.23	41.84	66.21	188.39	57.5	529.29	123.49

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₂.

Figure 2. Weighted mean depth (m) of Rotifera, Copepoda, and Cladocera biomass relative to the cline based on temperature and dissolved oxygen for the sampled subarctic thaw ponds

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₂ 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.

Figure 3. Vertical distribution of zooplankton biomass and the associated algal (phytoplankton biomass and EFA concentration) and nonalgal (Chaoborus and bacteria biomass) variables in KWK6, KWK12, and SAS2A ponds. Standard error bars are shown for variables with replicates. The horizontal dashed line defines the cline based on temperature and dissolved oxygen, between the mixed and stratified layers

In these three ponds, phytoplankton biomass varied from 0.1 to 70.2 mg L⁻¹ 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₂ concentrations, as was reflected in the model fit values of R² = 0.29 for total zooplankton, R² = 0.25 for rotifers, and R² = 0.16 for copepods (Table 3). Iterating Chaoborus, phytoplankton, and EFA into the models improved the fit, at R² = 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₂ alone (depicted by the AIC_c weight of 0.25), whereas the second (Chaoborus) and third (O₂ + 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₂, phytoplankton, and O₂ + phytoplankton; however in these models, the O₂-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². The other variables had a much smaller effect. For the copepods, O₂ 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² when integrated (Table 3). When considered alone, only O₂ was significant (DistLM marginal test, p < .05), explaining 29 percent of the overall variation of total zooplankton (Supporting Information Table S2).

Table 3. Best models explaining total zooplankton, rotifer, and copepod biomass distribution according to the distance‐based linear modeling (DistLM)

AICc	AICc weight	R^2	RSS	No. variables	Variable selection
Total zooplankton
78.09	0.25	0.29	3,539.4	1	DO
79.73	0.11	0.19	4,015.5	1	CB
80.12	0.09	0.36	3,170.5	2	DO, PB
Rotifera
81.17	0.31	0.25	4,484.8	1	DO
83.56	0.09	0.10	5,390.6	1	PB
84.09	0.07	0.28	4,302.3	2	DO, PB
Copepoda
91.90	0.17	0.16	10,238	1	DO
92.00	0.16	0.15	10,318	1	CB
93.21	0.09	0.29	8,678.1	2	CB, EFA

Note. Variables included in the models: CB = Chaoborus biomass, DO = dissolved oxygen, EFA = essential fatty acids, PB = phytoplankton biomass, RSS = residual sums of squares.

Discussion

We found large variation in temperature and O₂ 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₂ 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₂ 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₂, 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 (Rautio et al. 2011) 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 well-oxygenated water, the thaw pond zooplankton community composition favored a high abundance of rotifers. These organisms are known to be tolerant of low O₂ 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₂ 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 (Rautio et al. 2011), 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₂ 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₂ 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₂ 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₂ 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₂, 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₂ 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₂ 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⁻¹. 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₂ 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.

Acknowledgments

We offer our gratitude to the Cree and Inuit communities in Whapmagoostui–Kuujjuarapik for giving us access to their ancestral lands. We are also grateful to the Centre for Northern Studies (CEN) Whapmagoostui–Kuujjuarapik Station and its former manager, Claude Tremblay, for providing us with transportation and lab facilities. We thank Laura Forsström for her help and cheerfulness in the field as well as for the estimates of phytoplankton biovolume, Pierre Carrier-Corbeil for fatty acid analyses, Ginette Méthot for her expertise in zooplankton identification, and Murray Hay for his editorial assistance. Special thanks as well to Beatrix Beisner, Irena Creed, Jean-François Lapierre, and the anonymous reviewers for their judicious comments that improved this article.

Disclosure statement

No potential conflict of interest was reported by the authors.

Supplemental material

Supplemental data for this article can be accessed on the publisher’s website.

Additional information

Funding

Funding was provided by the Canada Research Chairs Program (Grant 950-217323), the Natural Sciences and Engineering Research Council of Canada, the Canada Foundation for Innovation, the Network of Centres of Excellence ArcticNet and the Centre for Northern Studies (CEN). The Ph.D. grant of Maxime Wauthy was also partly supported by the Merit Scholarship Program for Foreign Students from the Ministère de l’Éducation et de l’Enseignement Supérieur du Québec (Grant 193772).