The latest assessment report of the Intergovernmental Panel on Climate Change 2013 (IPCC AR5) summarized evidence for the warming of arctic and subarctic regions. Studies of the thermal state of permafrost in the southern Hudson Bay lowlands document air temperature increases of up to 2°C within the last 20 years (Smith et al., 2010). The date of sea ice break-up/freeze-up in Hudson Bay changed accordingly, leading to longer periods of open waters (Gagnon and Gough, 2005b) and large-scale environmental changes in areas along ecological transition zones (e.g., northern displacement of ecozones and tree-line) as assessed by satellite imagery time-series data (Piwowar and Ledrew, 1995; Fraser et al., 2012; Mamet and Kershaw, 2012). The densification of vegetation is another consequence of the rapid warming at high latitudes (Elmendorf et al., 2012). Such profound changes in the northern landscapes will alter ecosystem services and wildlife habitat resulting in modifications of species distribution and phenology, as well as biodiversity losses (Dokulil, 2014).
Lakes and ponds are numerous and widespread components of the North that respond sensitively to multiple environmental changes (Vincent and Pienitz, 1996). The limnological and hydrological conditions of surface freshwater ecosystems are known to change significantly across several marked ecozones within the subarctic region (SWIPA 2011) (Serreze et al., 2000). Detailed ecological surveys are therefore necessary to quantify changes in aquatic biota and their biodiversity, yet studies on aquatic ecosystems spanning multiple years are still rare (Lougheed et al., 2011). The Churchill Northern Studies Centre (CNSC, Churchill, Manitoba) is one of the few research facilities offering accessibility to the subarctic treeline and permafrost ecozone enabling the monitoring of ecosystems over longer periods of time. Natural sediment archives for reconstructing past environmental conditions over time can be tapped with paleolimnology (Smol, 2008). However, this approach often lacks the necessary temporal and spatial resolution to precisely unravel recent ecosystem changes (Prowse et al., 2006; Bouchard et al., 2013).
In order to document recent limnological and faunal changes in water bodies located across the boreal forest - tundra ecozone, we sampled and re-sampled a selection of lakes and ponds in the northwestern Hudson Bay Lowlands (HBL) in the vicinity of the CNSC in the summers of 1997 and 2006. This allowed us to compare the evolution of limnological conditions within the context of accelerated environmental change (Mamet and Kershaw, 2011, 2012). We also consulted historical ice maps, temperature and precipitation records to infer and quantify local environmental shifts over the period from 1978 to 2009. This information served to assess variations that occurred in the water quality and faunal assemblages between our two field surveys.
Freshwater microcrustaceans, in particular Cladocera and, to a lesser degree Ostracoda, are the most abundant species groups in Nordic water bodies and have been extensively used for water quality studies (Rautio, 2001; Sweetman and Smol, 2006; Wetterich et al., 2008; Rautio et al., 2011; Schneider et al., 2016). They are known to respond sensitively to environmental change and show distinct biogeographical distribution patterns (Jeffery et al., 2011; Rautio et al., 2011; Curry et al., 2012; Viehberg and Mesquita- Joanes, 2012). Previous biological studies on cladocerans and ostracods in our study area (Shelford and Twomey, 1941; Hebert and Hann, 1986; Havel et al., 1990a, 1990b; Boileau et al., 1992; Hann, 1995; Ng et al., 2009; Weider et al., 2010; Jeffery et al., 2011; Abbate and Sagri, 2012) enabled us to use this baseline information as a reference to compare against our most recent limnological survey completed in 2006. Our study aims to provide answers to fundamental questions related to modifications in faunal community composition, range (migration/dispersal) and biodiversity that are forced upon aquatic biota by the observed rapid ecoclimatic and environmental changes in the North.
Study region and sites
Air temperature and ice concentration data
The study region is located on the southwestern shore of Hudson Bay (N58.75° - 58.80°; W94.20° - 93.60°; Fig. 1) and is characterized by the transition of the northern limit of the boreal forest biome to the southern boundary of the arctic tundra biome. Churchill is one of the few cities in northern Canada that maintains a weather station since 1929 due to its (former) economic importance. The modern climate conditions in Churchill are characterized by four months with above 0°C mean air temperature (June, July, August and September). During this summer period about half of the mean annual precipitation (435 mm) is recorded (Environment Canada, 2013). The mean July air temperature varies between +8°C and +16°C, while the mean January air temperature varies between -34°C and -19°C. The mean annual air temperature averages -7°C. Consequently, lakes and ponds are covered with ice during eight months.
For our study, we used instrumental temperature data with monthly resolution for the last 32 years (1978 to 2009) from Churchill, Manitoba (Climate ID: 5060600 and 5060601, Environment Canada, 2013). For each year, we calculated the mean air temperature of the three months before ice freeze-up (August, September and October; ASO) and ice break-up (May, June and July; MJJ) as the length of this period is determined by ice conditions and consequently controls the duration of the life cycle of microcrustaceans (Fig. 2). For this purpose, weekly sea ice charts from Hudson Bay were re-evaluated for ice concentration during the same time period to estimate the yearly dates of ice freeze-up and ice break-up in SW Hudson Bay. To this end, we followed the methods described by Stirling et al. (1999) and Gagnon and Gough (2005b) and used a grid of 0.5° intervals south of 63° latitude and west of 90° longitude (91 grid points). The break-up of the annual sea ice was defined as the date by which half the grid points had a total ice cover of 5/10 or less, whereas the freeze-up was defined as the date by which half the grid points were equal to or exceeded 5/10 ice concentration (Etkin, 1991). Data on lake and pond ice cover are sparse and archived since 2003 by the staff of the CNSC (LA Fishback, personal communication). Past dates on lake or pond ice break-up can also be inferred from aerial photographs (A24922 #48-75; Department of Energy Mines and Resources - Natural Resources Canada, 1978) or published reports (McClure, 1943). Based on these results, we conclude that ice cover from tundra ponds generally breaks up 6 to 21 days earlier than on deeper and larger lakes in the same area, and that they eventually freeze up 2 to 9 days earlier. The ice cover of freshwater ecosystems follows the ice conditions on Hudson Bay by a lead of up to 48 days in the spring and up to 14 days in the fall.
The temperature and sea ice data provide evidence for increasingly warm mean temperatures in MJJ (+ 2.1°C) and ASO (+ 2.6°C) over the past 32 years, as well as periods of open water conditions on SW Hudson Bay that lengthened by 39 days from 122 to 161 days. Thus, we infer that tundra pond and lake catchments were also subjected to a similar extension in the length of the vegetation (growing) period.
Tundra landscape features
Abundant peatlands in open woodlands characterize the transitional landscape in the Hudson Bay Lowland (HBL). The peatlands are usually fens characterized by wet, dense, amorphous sedge and moss peats. Forested peatlands with thermokarst ponds can be found only on kame and esker sands, which are covered mainly by Black Spruce (Picea mariana) and Larch (Larix laricina) stands.
The surficial geology is made up of marine and glaciomarine sands and gravel deposits of the postglacial Tyrrell Sea (Dredge, 1992b). The only exception are local kame and esker sands that form prominent ridges of Late Wisconsinan age, as well as prominent ridges composed of aphanitic greywacke of Pre-Quaternary age along the coastline or Silurian dolomitic limestone further inland (Dredge, 1992b). Standing waters of different size are numerous and scattered throughout HBL and characterized by their shallow (<2.5 m depth) and alkaline (pH >8.5) waters, which can be explained by the underlying (glacio-) marine sediments and bedrocks consisting of Silurian dolomitic limestones (Dredge, 1992a).
The hydrological setting in HBL changes significantly throughout the year and was extensively investigated by Macrae et al. (2004) and White et al. (2014). During winter, the aquifers are inactive due to frozen soils and hence there is no connection between surficial waters. As air temperature (May) exceeds 0°C in the spring, the soil and aquifer layers are still frozen, but melt water forms nonchanneled overland flows that connect ponds and lakes. In the summer, the water levels of ponds and lakes are mainly controlled by precipitation and evaporative processes that can result in water-level fluctuations of up to 20 cm during periods of drought (Macrae et al., 2004). Strong precipitation events, however, may result in runoff and hydrological connection (White et al., 2014).
In total, 16 lakes and ponds were initially investigated during 29 July to 5 August 1997 by Claude Duguay (University of Waterloo, Ontario) and RP in the surroundings of Churchill, Manitoba, and across the transition zone of the northern treeline (Fig. 1). Nine years later, Ghislain Côté (Université Laval, Quebéc) and FV revisited eleven localities between 16 and 21 July 2006.
The tree-vegetation cover in the vicinity of each locality was classified according to three types, defined mainly by the presence of Picea mariana and Larix laricina along the lake shores. Shrub tundra-type (ST) is characterized by the absence of any tree vegetation or singular appearance of a tree along the shore. Forest-tundratype (FT; open woodland) was described when the tree stands were scattered and open area were occupied by shrubs and lichens. Boreal forest-type (BF; taiga) was defined when the locality was completely surrounded by closed tree stands.
The water samples were collected at the site of maximum water depth of each site. Profiles of water temperature, pH and specific conductivity were recorded each time with a multisensor Quanta Hydrolab®. Unfortunately, the oxygen sensor failed during the fieldwork campaign in 2006, but first readings presume that all sites were well oxygenated. For chemical analysis, water samples were taken with a Kemmerer sampler and stored in PE-bottles in a cooler box until the samples were split and pre-treated in the laboratories of the CNSC within 12 hours according to the protocol of the Analytical Methods Manual of Environment Canada laboratories. The samples were then stored at 4°C, shipped to and analyzed by the National Water Research Center (Burlington, Ontario) successfully for 10 parameters: chlorophyll-a (Chla), dissolved organic carbon (DOC), dissolved inorganic carbon (DIC), total Kjeldahl nitrogen (TKN), total phosphorus (TP; filtered water), sulfate (SO4), chloride (Cl), silicate ions (SiO2), calcium (Ca), and sodium (Na).
Zooplankton samples were retrieved quantitatively in 2006 by three vertical tows of a conical plankton net (63 μm mesh size) with a mouth diameter of 40 cm from 1.0 m water depth to surface close to the center of the sampling site. In shallow waters, six vertical hauls of 50 cm were sampled instead. The total water volume sampled encompasses approx. 380 L. The captured material was fixed in the field with 70% ethanol and stored cool.
A set of seven surface sediment samples was collected with an Aqua Research® gravity corer with a discrete sampling area of 44.2 cm² (total 309.4 cm²). The material was used to retrieve abundance data from meiofauna (i.e., ostracods and benthic cladocerans) and organic sediment analyses at each site. It is noteworthy that in some ponds, gravel, stones and boulders dominated the sediment and thus the sampling position was selected based on the presence of penetrable sediment patches. The supernatant water and the top 2-3 cm of sediment from the short cores were conserved with 70% ethanol and stored cool. For each site, the samples were taken from different positions to ensure a representative minimum area for the coverage of the species diversity (Viehberg, 2006).
In the Aquatic Paleoecology Laboratory (APL) at Université Laval, we used a Folsom splitter (McEwen et al., 1954; Van Guelpen et al., 1982) for zooplankton samples when the abundance exceeded 500 individuals per sample. The results of the zooplankton counting were standardized to a water volume of 100 L. Sediment samples were washed over a 125 μm sieve gently under warm water. Ostracods and cladocerans were isolated from the sediment samples by using a flotation technique and an exhauster device described by Viehberg (2002). All sediment samples were post-screened for remaining individuals using a stereo binocular. Individuals were counted as living species when the soft-body was attached to their carapaces. The abundance of benthic microcrustaceans was standardized to 100 cm². The identification of cladocerans and ostracods followed Benzie (2005), Brooks (1959), Delorme (1969, 1970a, 1970b, 1970c, 1970d, 1971, 2001), Frey (1959, 1962), Meisch (2000), Orlova- Bienkowskaja (1998) and Tressler (1959).
For sediment characterization, a subsample of 15.6 to 25.2 g of wet sediment was analyzed for water content and loss-on-ignition (LOI), following the instructions outlined by Heiri et al. (2001). The sediment was weighed before and after drying the sediment at 105°C in the oven overnight to assess the water content. Then the organic matter was estimated by weighing the loss of dry sediment after burning the samples at 550°C for 8h in crucibles. Finally, the inorganic fraction was assessed by burning the previous samples at 950°C for another 8 h and monitoring the total loss.
Numerical analysis and diversity indices
All statistical tests and unconstrained ordination techniques were completed with the software package R version 3.2.5 including the libraries: ade4, vegan (R Core Team, 2016). A paired sample t-test was used to identify significant changes (P<0.01) between environmental variables measured at identical sample sites during the campaigns of 1997 and 2006. In addition, we were interested in the consistency of the total variance of each environmental variable between the two fieldwork campaigns, therefore we used an one-way ANOVA approach to test on significance. A principal component analysis (PCA) was applied to standardized environmental variables to identify driving factors characterizing the ponds and lakes in our study area and to assess the multivariable change over time. To explore the main spatial and ecological dimensions of species distribution and studied sites, we ran a correspondence analysis (CA). The results are scaled by the weighted average of species abundances (scaling 2) to focus on the ordination of species. The influencing environmental variables whose eigenvalues contribute more than average to the variation in our data set were projected posteriori in the CA biplot by regression on the two axes. A permutation test (n=999) assessed the significance of the regression coefficient (r²). To express the species diversity of the studied lakes and ponds (alpha diversity), we calculated for each site the species richness as the sum of all different species occurrences and converted the entropy to diversity measures (Jost, 2006). The Shannon entropy index (1) was calculated as follows:
where pi is the relative abundance of the ith species focusing on the abundant species. For direct comparison reasons, we converted the entropy index (1) to a diversity measure, the Shannon diversity index
and similarly, we calculated the reciprocal of Simpson diversity indices as
to evaluate the diversity in terms of dominating species (Hill, 1973; Jost, 2006). The indices were determined for zooplankton and benthos separately, because the results of the different sampling methods are incomparable. To evaluate the variation of species assemblages between the lakes and ponds, we used two measures, the Sørensen (βsor) similarity index:
and the Simpson (βsim) similarity index:
where a is the number of species both sites have in common, b and c are the number of species that are exclusive in one or the other (sensu Lennon et al., 2001). Both indices have a lower limit of zero when the sites have no species in common and an upper limit of one when the sites have an identical species inventory. The βsor is more sensitive to detect total differences in richness between sites, while βsim better reveals compositional differences even if the richness is rather asymmetric between sites.
RESULTS AND DISCUSSION
The shallow water basins of our study lakes and ponds are representative of tundra lakes from throughout the circumpolar subarctic landscapes; they were well oxygenized at the bottom and not stratified (Smol et al., 2005) (Fig. 3). The analyzed chemical parameters indicate that all investigated lakes are oligotrophic and dilute (Tab. 1). However, due to their limited volume, the physical and chemical variables are directly influenced by climate forcing mechanisms (e.g., precipitation and temperature) and the development of their catchment (e.g., geochemical leaching from the bedrock, local vegetation, peat accumulation and erosion) (Macrae et al., 2004; Pienitz et al., 2008). In particular, the water temperature is coherent with the air-temperature regime as well as solar radiation. On a daily scale, temperature changes of up to 4°C were observed (Fig. 3; profiles). On a yearly scale, the date of ice break-up on ponds and lakes in the study is respectively 31 to 48 and 27 to 45 days earlier than that of the sea ice in western Hudson Bay (LA Fishback, personal communication; Fig. 2). It has been shown previously that the distinct dates are strongly related to the mean air temperature of the months May, June and July (r²=0.33; P<0.001; this study) (Gough et al., 2004; Gagnon and Gough, 2005a, 2005b). Two additional references to ice break-up dates of freshwater ponds in the Churchill area are reported by McClure (1943) as day 162 in 1937, whereas a set of aerial images (National Air Photo Library, Natural Resources Canada, Film A24922) from 22 June 1978 documents the approximate date of ice breakup on lakes (day 173) in our study area. This latter data point corresponds well to the historical temperature record as shown in Fig. 2.
A significant positive correlation (r²=0.48; P<0.001) also exists between mean air temperatures of August, September and October and the date of ice freeze-up (Fig. 2). The time of ice freeze-up of the ponds and lakes is simultaneous or up to 33 to 48 days earlier than the formation of sea ice in western Hudson Bay. As a result, the time of open water and increased biological production varies between 118 and 165 days in ponds and 100 and 156 days in lakes, thereby influencing also the development of the microcrustacean fauna (Wetterich et al., 2008). In this study, we conclude from the measured data that the average monthly winter air temperature increased by approx. 2°C during the investigated time span of 30 years (Fig. 2), while average monthly precipitation over the past decades increased by approx. 5 mm (Environment Canada, 2013). In case this trend is continuous, we anticipate an increase of ice-free waters in ponds and lakes by approx. 12 days per decade.
Historical changes in water quality
The general hydrology of ponds and lakes in the study area varies substantially after pronounced precipitation and/or snowmelt events (Macrae et al., 2004). To compare data, it is important to highlight the meteorological setting before and during the fieldwork campaigns in 1997 and in 2006. The mean air temperature during the campaigns was comparable, 15.5°C (min. 8.0°C; max. 28.5°C) in 1997 and 15.6°C (min. 8.8°C; max. 24.0°C) in 2006. In contrast, the precipitation that accumulated during the week preceding the start of the fieldwork missions in 1997 and 2006 was different. Precipitation totalled 19.5 mm and 82.5 mm in 1997 and 2006, respectively, but no major rain event occurred during both sampling periods. The measured specific conductivity of the ponds was systematically and significantly (P=0.004) lower in 2006 than in 1997, but the general trend shows increased values with proximity to the Hudson Bay shoreline (Tab. 1), likely caused by sea spray inputs (Pienitz et al., 1997). As a result, the ponds within a distance of 2 km have an increased specific conductivity between 252 and 563 μS cm–2, while distant ponds are more dilute between 100 and 148 μS cm–2. As a consequence, the main variance in the PCA is associated with the conductivity (salinity) gradient (PCaxis 1; eigenvalue=6.54; explained proportion=54.5%) and showing the dilution effect in water bodies revisited in 2006 (Figs. 4 and 5).
The PCA-analysis also reveals changes that evolve along PC-axis 2 (eigenvalue=1.83; explained proportion= 15.2%) (Fig. 4) that is loaded mainly by pH variations (score=0.94). Almost all studied water bodies are affected by this trend, with the exception of sites 9 and 12. Correspondingly, we measured systematically and significantly higher pH values in pond waters in 1997 (mean 8.7, sd. 0.2; Fig. 5) compared to field data in 2006 (mean 8.0, sd. 0.2; Fig. 5). We believe that the resulting lower pH values are due to several causes. Heavy precipitation events prior to the 2006 field campaign must have led to acidification and general ion-dilution in our studied waters as mentioned above (Galloway et al., 1984).
A slight eutrophication of the pond waters over the years may be interpreted based on the higher mean DOC and Chla concentrations in 2006, while nutrients and silicate appear to have been metabolized and taken up by aquatic biota at the time of sampling (Fig. 5). Other studies have shown the fertilizing impact of DOC on primary productivity in northern Swedish lakes (Seekell et al., 2015). Future investigations may test this hypothesis by extending the seasonal coverage of limnological surveys and by putting special emphasis on the impact of prolonged ice-free periods on lake trophic changes.
Species inventory - Alpha diversity
In our samples from 2006, microcrustaceans are the dominant meiofauna group. For reasons of comparison, we focused on Cladocera and Ostracoda as both groups were subject of previously published (partly molecular biology) studies in HBL (Shelford and Twomey, 1941; Hebert and Hann, 1986; Havel et al., 1990a, 1990b; Boileau et al., 1992; Hann, 1995; Weider et al., 1999a, 1999b; Marková et al., 2007; Ng et al., 2009; Xu et al., 2009; Weider et al., 2010; Jeffery et al., 2011). In addition, both microcrustacean groups have potential value for paleoenvironmental assessment studies and this study is intentionally necessary to interpret future analyses of faunal assemblages in the HBL.
We identified a total of 44 species (Tab. 2; Cladocera: 18, Ostracoda: 26), of which nine species are known to be planktic (sensu strictu) (Dodson and Frey, 2001). Our sampling effort resulted in eleven species that were previously not reported from HBL waters (Alona (A.) circumfimbriata, A. guttata, A. quadrangularis, Candona (Can.) acuta, Can. decora, Cypria ophtalmica, Graptoleberis sp., Limnocythere (L.) inopinata, L. itasca, L. liporeticulata and L. varia). On average, our study sites comprised 21 species (Benthic: 14, Planktic: 7). Due to the shallow water depth, we collected several benthic chydorids in the plankton net samples. However, the abundance data from the parallel sediment samples were less variable. Therefore, they are more reliable and representative to interpret species occurrences and to calculate diversity and similarity indices (Tabs. 3 and 4). It is noteworthy that the species richness is in favour of the benthic species assemblage, mainly because the complete ostracod fauna is exclusively benthic and comprises 26 species.
The most abundant and common planktic species was Daphnia (D.) tenebrosa (DATE), which occurred in all studied sites and in high abundance particularly in waters with higher concentrations of Chla (i.e., high phytoplankton concentration) of up to 950 n/100L (Fig. 6). Similarly, other abundant planktic species, D. middendorffiana (DAMI) and Holopedium gibberum (HOGI), also show a positive correlation to increased concentrations of Chla (Fig. 6). In three sites (4, 8 and 14), all planktic species were found, but their faunal assemblage was dominated by D. tenebrosa and D. middendorffiana thus resulting in relatively low diversity indices (Tab. 3). In contrast, a minimum of only three planktic species were found in sites 9 and 12 and thus resulting in the minimal diversity indices of D1Plankton=1.11 or 1/D2Plankton=1.04 and D1Plankton=1.49 or 1/D2Plankton=1.29, respectively. An explanation for the low diversity in planktic Cladocera could be limited food availability inferred from low Chla concentrations and increased DOC and TP concentrations (Tab. 1). Another explanation may be due to competition by Copepoda (Gliwicz and Umana, 1994), which were not in focus of this study. Highest diversity indices were found in sites located within the boreal forest zone (6, 7, 13 and 14) ranging between D1Plankton=2.88-4.93 and 1/D2Plankton=2.26-4.05, respectively.
Alonella excisa (ALEX) and Alona guttata (ALGU) are the most common benthic microcrustaceans, which were caught in all studied sites with abundances of up to 56 n/100cm². Followed by Picripleuroxus striatus, Chydorus sphaericus, Cyclocypris (Cyc.) globosa, Cyc. ovum and Pseudocandona (Psc.) compressa (Tab. 2) that are found in the majority of the sites. Again, site 9 holds the lowest species richness with a total of nine benthic species and thus resulting in overall low diversity indices in benthic microcrustaceans: respectively D1Benthos =4.52 or 1/D2Benthos=2.97. The neighboring sites and measured variables give no reasonable explanation for the minimal diversity. In general, sites in the boreal forest zone (6, 7, 13 and 14) yield a low number of benthic species (Tab. 3), whereas waters with higher conductivity and closer proximity to the coastline tend to have higher diversities up to D1Benthos=14.09 or 1/D2Benthos=12.36 (Tabs. 2 and 3), while site 9 is an exception.
Axis 1 of the correspondence analysis (CA) explains 34.0% of the total inertia (0.660) in the dataset and this axis is highly loaded with inland sites and their species inventory, especially ostracod species, such as Bradleystrandesia reticulata (BRRE), Can. candida (CACA), Can. decora (CADE), Can. distincta (CADI), Cyc. sharpie (CYSH), Limnocythere (L.) itasca (LIIT) and Psc. rostrata (PSRO) (Fig. 6). As expected, the projection of the environmental variable distance from the coast (DIST) is significantly correlated (r²=0.71; P=0.002) to the CA results.
CA Axis 2 explains a further 19.5% and is loaded with sites that are characterized by their high Chla concentrations (r²=0.53; P=0.026) and/or high LOI550 values (r²=0.43; P=0.034). In particular, the abundance of Holopedium gibberum (HOGI) is positively correlated to higher Chla concentrations. In contrast, benthic species such as Cyp. ophtalmica (CYOP), Gryptoleberis sp. (GRSP), Cyc. ampla (CYAM), Eurycercus glacialis (EUGL) and Bosmina sp. (BOSP) seem to favour sites with organic-rich sediment type (increased LOI550).
Site evolution – Beta diversity
The similarity matrix enabled us to compare the species inventory between sites more generally (Tab. 4). In the study area, sites in close geographical vicinity often show high βsor or βsim values such as neighboring sites 6 and 7 (βsor=0.703, βsim=0.765) or the cluster of sites 2, 3 and 4 (βsor=0.837-0.636, βsim=0.947-0.737). Lower values found in the cluster of sites 8, 9 and 10 (βsor=0.720-0.410, βsim=0.917-0.667) are partly caused by the overall low diversity in site 9, as discussed above.
Focusing on the closest similarities, sites located within the boreal forest (6,7,13 and 14) show a distinct cluster (βsor=0.727-0.524 and βsim=0.800-0.588) while the fauna of site 13 is least similar in this context, partly because of the overall low species diversity of benthic microcrustaceans at this site.
A further cluster of sites with similar species assemblages are sites close to the HBL coastline (2, 3, 4, 8, 10) with βsor=0.837-0.618 and βsim=0.947-0.684. It is noteworthy that these sites also show high alpha diversity indices (Tab. 3). According to previous reports on microcrustaceans from HBL, the species richness in the coastal rock pools is higher than in the present study of lakes and ponds, including 34 species that were not sampled in the present study (Havel et al., 1990a,1990b; Jeffery et al., 2011). We suggest two main causes for the substantial differences in species richness: i) The dramatic variation in salinity 150-13,650 μS/cm–1 in rock pools expresses a greater ecological gradient to host a more diverse microcrustacean fauna (Weider and Hebert, 1987; Havel et al., 1990a, 1990b; Weider et al., 2010); ii) The intense sampling effort of numerous rock pools by various studies focussing on cladocerans covers a time span of 25 years and includes also rare species (Weider et al., 2010). Therefore, a direct comparison of species diversity or similarities is a matter of scaling (Chase and Knight, 2013).
Sites in boreal forest – gamma diversity
The knowledge of the (genetic) biodiversity and ecology of Arctic and Subarctic Cladocera and Ostracoda species is based on numerous studies that were motivated by different research goals, starting with biological expeditions near the end of the 19th century (Sars, 1899). We present a species list of studies completed in the Subarctic and Arctic realms of North America (Supplementary Tabs. 1 and 2). Surprisingly, we found no evidence that both Cladocera and Ostracoda were investigated in the same study. Hence, to assess the species diversity of the studied lakes and ponds, we compared the fauna with two surveys which investigated separately the ecology of cladocerans and ostracods in lakes and ponds in the boreal forest biome in southwest Yukon with similar hydrological settings in terms of the physical and chemical water quality parameters (Swadling et al., 2000; Bunbury and Gajewski, 2005). The species richness in total numbers is comparable (Tab. 3), while the relatively low similarity indices (βsor=0.420 and βsim=0.459) indicate distinct biogeographical distribution of individual species (Tab. 2; Adamowicz et al., 2009; Curry et al., 2012). Weider and Hobæk (2003) suggested that the asynchronous deglaciation history of the Cordilleran and Laurentide ice sheets resulted in different recolonization patterns of the species in the Arctic and Subarctic that caused species dependent phylogeography. Thus, it appears that proglacial aquatic ecosystems in the Yukon evolved when the landscape became ice-free already at ca. 13,000 years BP (Dyke and Prest, 1987; Bunbury and Gajewski, 2009b; Bunbury, 2012). In contrast, HBL deglaciated and was submerged by marine waters of the postglacial Tyrrell Sea around 8000 years BP. However, freshwater ecosystems established not earlier than 5000 years BP as the landscape gradually emerged above sea-level through glacio-isostatic rebound (Dyke and Prest, 1987; Dredge, 1992b).
We showed that lake and pond ecosystems in HBL are exposed to similar climate stresses as shown also for arctic tundra ponds (Smol et al., 2005). The waters including the rock pools are identified to be vulnerable, as they have little to no buffer capacity to changes in temperature and precipitation rates (Weider et al., 2010; White et al., 2014). The evident climate warming will shorten the period of ice cover and inversely extend the vegetation period, which will impact phyto- and zooplankton life cycles, but also alter the hydrology due to increased evapotranspiration and absorption of solar radiation (Arp et al., 2015). Consequently, we will experience the densification of the vegetation in the catchment, i.e. the northward displacement of the northern treeline and dispersal of boreal (aquatic and terrestrial) species and a retreat of arctic and subarctic (aquatic and terrestrial) species (Jacques et al., 2016).
In this context, long-term ecological research (LTER) is essential to detect systematic floral or faunal changes and assess the resilience, reactivity, and variance of ecosystems based on abiotic and biotic factors (Haase et al., 2016). Our assessment of the faunal similarity and diversity of cladoceran and ostracod species in freshwaters within the boreal forest and shrub tundra catchments revealed distinct faunal assemblages. Markedly, only certain ostracod species (i.e., Can. acuta, Can. acutula, Can. decora, Can. distincta and L. liporeticulata) appear exclusively in boreal forest sites. Considering also previous studies (Supplementary Tabs. 1 and 2), we suggest especially Can. acuta, Can. acutula and Can. decora as indicator species of proceeding warming and eventually following the northern treeline and identified Can. rawsoni and Ton. glacialis as subarctic and arctic faunal elements. Both species were reported also from other Arctic sites (Røen, 1962, 1968; Wojtasik, 2008; Bunbury and Gajewski, 2009a). Despite single indicative species, overall it appears that waters in boreal settings have a higher alpha diversity in zooplankton, whereas waters in the shrub tundra display a higher diversity in benthos faunal assemblage. Further research is needed to monitor these trends and to specify implied consequences to the food web (Thienpont et al., 2015).
We conclude that our study contributes quantitative results that are important to understand the decadal hydrochemical and species variation of subarctic waters that are extremely abundant in the northern Manitoba landscape. In combination with other regional reference datasets (Jacques et al., 2016), the resilience, reactivity, and variance of these aquatic ecosystems may be assessed in detail to predict ecosystem changes. Meanwhile, complementary palaeolimnological studies (including remains of diatoms, cladocerans and ostracods) may hindcast the ecosystem dynamics of the past.