Northern regions are expected to experience large environmental change over the next few decades. The response of biota will depend on changes in the local environment, regional processes that influence lake connectivity, and species interactions. In 2008, we surveyed 92 lakes and ponds across Wapusk National Park, located on the southwestern shore of Hudson Bay. At each site we assessed water chemistry and zooplankton community composition. In an effort to understand how the aquatic ecosystems will respond to future environmental change, we determined local characteristics (e.g., water chemistry), regional spatial factors (e.g., dispersal), and biotic interactions (e.g., species associations) influencing community composition. Important environmental variables included lake area, pH, ionic composition, total phosphorus, and chlorophyll a; however, spatial variables explained more variation than environmental variables, suggesting that dispersal is an important driver of zooplankton composition in this region. Additionally, species exhibited negative co-occurrence patterns, suggesting biotic interactions are important in structuring the zooplankton communities. As environmental conditions change and the distribution of habitat (i.e., coastal fen, interior peatland, and spruce forest) shifts, evidence that the park's zooplankton community is spatially structured coupled with our suspicion that zooplankton are likely to experience high dispersal levels in Wapusk leads us to suggest zooplankton may indeed be able to track changing environmental conditions within the park, although it remains unclear how species interactions will modify this expectation.
One of the primary goals of community ecology has long been the identification of the factors that control the distribution of species (Hubbell, 2001), a task imbued with particular importance given predictions of future environmental change (Meehl et al., 2007). Species distributions, and the ecosystem functions that species provide, are likely controlled by environmental factors that determine which ecological niches are available, and spatial factors that determine species' access to those niches (Leibold et al., 2004). Furthermore, both classes of factors are expected to play a role in determining how ecosystems respond to environmental change, as the extent to which ecosystem functions will be affected by environmental change likely depends on the ability of a biota to adapt to new environmental conditions (Bell and Gonzalez, 2009), as well as its ability to disperse and thus track those environmental conditions (Leibold and Norberg, 2004; Bell and Gonzalez, 2011).
Wapusk National Park (WNP), located on the southwestern shore of Hudson Bay near Churchill, Manitoba, provides an attractive system to study this interaction between environmental conditions and dispersal capabilities of aquatic organisms. Water is a prominent feature on the landscape, covering ∼40% of the 11,475 km2 total area of WNP. Within the park, there are over 10,000 lakes/ponds embedded in regions differing in lake density and connectivity. Canada's subarctic is expected to experience dramatic changes in temperature in the near future (Christensen et al., 2007) and is generally a region where biota are poorly understood and in need of long-term monitoring (Rouse et al., 1997). The park itself contains a wide range of environmental conditions (Symons et al., 2012), which creates strong potential for niche differentiation among different species, and the park's high concentration of water bodies provides the potential to examine the importance of spatial factors in a habitat likely to have high dispersal potential for aquatic organisms, such as zooplankton.
Zooplankton are a diverse group of aquatic organisms adapted to a wide variety of environmental conditions, and they play an important role in trophic energy flow through their ecosystems (Strecker and Arnott, 2008). It has been long established that zooplankton have strong dispersal capabilities, largely due to the potential for many members of the group to form stress-tolerant resting stages (e.g., Proctor, 1964). Consequently, it is widely hypothesized that zooplankton distributions are closely tied to environmental conditions (Leibold et al., 2004), and, presumably, the species' fundamental or realized niches (Hutchinson, 1957). Indeed, this link between environmental conditions and zooplankton distributions has some strong experimental and field support (e.g., Cottenie et al., 2003; Cottenie and De Meester, 2004; Beisner et al., 2006; Strecker et al., 2008; Gray et al., 2012).
Conversely, it has also been hypothesized that various factors are capable of disrupting the relationship between biota and niche optima. For example, high dispersal conditions may create masseffects situations where species disperse so rapidly that they are present even in suboptimal habitats (Leibold et al., 2004), and in habitats where communities are already well established, strong preexisting biotic interactions may prevent new species from establishing (Shurin, 2000). Indeed, the potential for such priority effects is high in zooplankton, given their rapid generation times and potential to form “banks” of resting stages (De Meester et al., 2002). Furthermore, the neutral perspective hypothesizes that some species distributions may not reflect differing responses to environmental variation, but rather that community composition may be a result of both stochastic competitive dynamics among roughly equivalent member species, and spatial factors that affect species access to habitats (Hubbell, 2001). There is a near consensus now that ecological equivalence is unlikely for broad groups of organisms; however, there is evidence that small groups of species may be ecologically similar, causing stochasticity to be a major factor driving community dynamics (Vellend, 2010).
The metacommunity paradigm now popular in ecology represents a synthesis position which recognizes that environmental conditions likely have some impact on species distributions, but that the influence of the environment is tempered by spatial influences, be they low dispersal that prevents species from accessing habitats to which they are well suited, or high dispersal that maintains species in habitats to which they are poorly adapted. Here we examine a large data set of zooplankton communities over a marked environmental gradient in Canada's subarctic that is well suited to analyze the extent to which zooplankton distribution is a result of environmental conditions, spatial structure, competition, or some mix of all of these factors. We had three goals: (1) to reveal important environmental drivers of zooplankton community composition; (2) to determine the relative importance of local environmental conditions compared to spatial variables (i.e., dispersal); and (3) to evaluate the extent to which species interactions shape local community composition. These analyses will help inform on the extent to which changing environmental factors are likely to affect zooplankton communities, and thereby influence broad-scale ecosystem function, in a region expected to undergo significant change.
WATER AND ZOOPLANKTON SAMPLING
From 22 July to 3 August 2008 we used a helicopter to sample 92 lakes/ponds in Wapusk National Park. These lakes/ponds were distributed throughout many different types of land cover, which have been grouped into three main habitat types: spruce forest, interior peatland, and coastal fen. Spruce forest is used to describe areas of lichen spruce bog and regenerating burn areas; interior peatland is used to describe areas of lichen melt pond bog and lichen peat plateau bog; and coastal fen describes areas of sedge bulrush-poor fen, sedge-rich fen, and sphagnum larch fen (land cover classification from Brook, 2005). Water bodies were named arbitrarily by the order they were sampled and latitude/longitude, and physical/chemical characteristics are listed in Appendix Tables A1–A3. Most lakes/ponds were small, having depths usually <0.5 m and surface area between 0.05 and 945 ha (median = 4 ha). Near shore at each site, we took in situ measurements of temperature and dissolved oxygen using a YSI 600OMS probe (YSI Incorporated, Yellow Springs, Ohio, USA). Weather was variable during the sampling period, and water temperature was correlated with the daily average air temperature (Pearson correlation, r = 0.21, p = 0.04). Water samples were collected 10–20 cm below the surface for water chemistry and chlorophyll a analyses. Water was filtered through a 75-µm mesh to remove zooplankton. For the analyses of major ions (Cl, SO4, Ca, Mg, K, Na, SiO2), nutrients (total phosphorus [TP], filtered and unfiltered, NO3-NO2, NH3, total nitrogen [TN]), alkalinity, dissolved organic carbon (DOC), dissolved inorganic carbon (DIC), specific conductivity, and pH, water samples were sent to Environment Canada's National Laboratory for Environmental Testing at the Canada Centre for Inland Waters (Burlington, Ontario) (EC, 1994). DIC, DOC, and filtered TP samples were filtered through a 0.45-µm celluloseacetate filter. TP samples were preserved with 1 mL of 30% H2SO4 prior to shipping. For chlorophyll a analysis, a known volume of water was filtered through a glass fiber filter (Whatman GF/C, pore size 1.2 µm), which was frozen until extraction. Chlorophyll a concentration was then determined using a Turner Designs TD700 fluorometer following a 24-h methanol extraction. Qualitative zooplankton samples were collected from each lake by horizontally towing a 75-µm mesh conical net through the water for approximately 10 m. Samples were immediately preserved in 70% ethanol for future identification.
Crustacean zooplankton were enumerated using a Leica MZ16 microscope at 40× magnification and were usually identified to species, except Alona spp. and chydorids, which were identified to genus; copepod juveniles and harpacticoids, which were identified to order; and ostracods, which were identified to subclass. A target of 200 individuals for each of crustaceans and rotifers, not including juveniles, was processed in successive subsamples of a defined volume. All crustaceans were identified using the keys of Brooks (1959), Wilson and Yeatman (1959), Smith and Fernando (1978), Dussart (1985), De Melo and Hebert (1994), and Hebert (1995). Rotifers were identified to genus for monogonont rotifers and to class for bdelloid rotifers. All rotifers were identified using Edmondson (1959) and Stemberger (1979). Given that we used 75-µm mesh to sample the zooplankton community, we are only considering the subset of the rotifer community >75 µm in size.
Lake/pond surface area and proximity to other lakes/ponds were assessed using ArcGIS 9.3. Lakes/ponds were derived from cloud-free Landsat 7 imagery from August 2001.
CANONICAL CORRESPONDENCE ANALYSIS
A multivariate ordination technique, canonical correspondence analyses (CCA), was used to determine if zooplankton taxonomic composition was related to environmental variables. First, unidentified monogononts and all insect/ arachnid/cnidarian taxa were removed. Second, extremely rare species were removed (>5% occurrence) as they can have a disproportionate effect on ordination results (Quinn and Keough, 2002). Exceptions were made for Bosmina freyi, which was lumped with Bosmina liederi to create the new taxon Bosmina spp., and for Ilyocryptus acutifrons, which was lumped with Ilyocryptus sordidus to make Ilyocryptus spp. Zooplankton relative abundances were Hellinger-transformed to reduce the influence of rare species, and zeros that are common in community data (Legendre and Gallagher, 2001). A CCA was appropriate because gradient lengths were long (>3) (ter Braak and Šmilauer, 1998). Correlations between environmental variables were investigated using variation inflation factors (VIFs). The ion-related data were highly correlated, so a principal component analysis (PCA) was completed using a correlation matrix of altitude, DIC, alkalinity, calcium, log(magnesium), log(specific conductivity), log(sulfate), log(silicate), log(potassium), log(sodium), and log(chloride). The first principal component of this PCA explained 90.3% of the variation in ion-related data between lakes/ponds and was used as an environmental variable for future analyses. After completing the PCA, all VIFs were <10, therefore all variables were retained (Quinn and Keough, 2002). The environmental variables, transformed where desirable to improve normality, included log(TN), log(TP), log(Area), log(Chl), log(DOC), temperature, pH, and the first principal component axis of the ion PCA. Unfortunately, we have limited information on predators in the lakes/ponds. Two of the largest lakes (Mary and Lee Lakes) have fish based on reports of fishing; however, the majority of the survey ponds are likely fishless because they freeze to the bottom in winter, owing to their shallow depth. Environmental variables were forward-selected using Monte Carlo permutation tests at p < 0.05 with 999 iterations. After variables were selected, Monte Carlo permutation tests were used to determine the significance of constrained axes. The DCA was completed using CANOCO 4.5 (ter Braak and Šmilauer, 2002), and the CCA was completed using R (R Development Core Team, 2012).
VARIATION PARTITIONING OF ENVIRONMENTAL AND SPATIAL VARIABLES
To assess the relative contribution of local environmental conditions and the spatial arrangement of lakes/ponds in the landscape to determining zooplankton community structure, we used a combination of spatial modeling and multivariate ordinations. Spatial modeling of study lakes/ponds was conducted using Moran's eigenvector maps as outlined by Dray et al. (2006). The spatial arrangement of lakes/ponds on the landscape was modeled using four different connection schemes to create connectivity matrices—Delauney triangulation, Gabriel graphs, relative neighbor graphs, and sphere of influence graphs (Dray et al., 2006). Each of these frameworks determines the proximity of lakes/ponds to each other and expresses these distances in a pairwise connectivity matrix. For each framework, the Hadamard product of each connectivity matrix by a variety of spatial weighting matrices (linear = 1 - dij/dmax, concave down = 1 - dija y/ dmaxa y, and concave up = 1/dijb y, for y = 1:20) was calculated to determine the potential spatial weighting matrices, and the model with the lowest Akaike information criterion (AIC) from this set was selected as the best. We then reran the AIC selection process to identify variables to be retained in the model after we had removed all eigenvectors not correlated with Moran's I, as these likely do not show spatial structure (Dray et al., 2006). The generation of the spatial predictors was performed using the R packages spdep (Bivand et al., 2012), spacemakeR (Dray, 2010), and their associated packages (R Development Core Team, 2012).
The environmental variables and zooplankton community composition data used in the ordinations for the variance partitioning analyses were the same as the CCA.
Variance partitioning was performed on the zooplankton community composition at different spatial scales within the park. As we felt that rotifers, cladocerans, and copepods were all likely to experience different spatial and environmental regimes, we performed three separate variance partitions, one for each group. Variance partitioning was performed in R using the varpart function of the package vegan in R (Oksanen et al., 2012; R Development Core Team, 2012). This function uses the species, environment, and spatial matrices for redundancy analyses (RDAs) and partial RDAs to calculate R2 adj values that represent the independent and shared variance explained by environment and space. A CCAbased variation partitioning method would be more optimal for our data due to long gradient lengths; however, there are no methods to calculate the adjusted R2 values for CCA-based variation partitioning (Oksanen et al., 2012).
ANALYSIS OF SPECIES CO-OCCURRENCE PATTERNS
To determine if biotic interactions are important in structuring local community composition, we completed a co-occurrence analysis that determines if there are negative (i.e., segregations) or positive (i.e., aggregations) species associations using the program EcoSim v. 7.0 (Gotelli and Entsminger, 2009). Species data were organized as a presence-absence matrix, where each row was a different species and each column was a different site. This analysis was completed for the total zooplankton community and the rotifers, cladocerans, and copepods separately. The C-score (Stone and Roberts, 1990) was used as an overall measure of species cooccurrence, as it quantifies the average number of checkerboard units (CU) between each species pair. A checkerboard is a submatrix of the form:
Therefore, a checkerboard unit represents an instance of negative co-occurrence between two species. To compare between the four different analyses (total zooplankton, rotifers, cladocerans, copepods), the C-score was standardized to account for differences in matrix dimensions. The standardized effect size (SES) = (observed index — mean of simulated indices)/standard deviation of simulated indices. Negative SES values represent positive species co-occurrences, or aggregations, and positive values represent negative species co-occurrences, or segregations (Stone and Roberts, 1990).
To determine the significance of the observed C-scores, Monte Carlo randomizations of community presence-absence data were used to create “pseudo-communities.” We used fixed-fixed constraints on the randomized matrices—i.e., row and column constraints that maintain row sums (number of occurrences of each species in the data set remains the same) and column sums (number of species in each site remains the same). This reduces the chance of Type I error (Gotelli and Entsminger, 2009) and is well suited to data recorded across heterogeneous environments (Rooney, 2008). The randomization occurred using an independent swap algorithm in which the original data matrix is shuffled by swapping random checkerboard submatrices after discarding the first 50,000 swaps. C-scores were calculated after 5000 randomized matrices were generated. The average C-score of the simulated matrices is included in the calculation of SES. All SES that are >1.96 or ≤1.96 are significant at p < 0.05 and suggest that negative or positive associations are different from potential random patterns (Gotelli and Entsminger, 2009).
DETECTING PROCESSES UNDERLYING SPECIES COOCCURRENCE PATTERNS
Ulrich and Gotelli (2013) showed that when there are multiple patterns in the matrices (i.e., both aggregations and segregations) the results of co-occurrence analyses can be misleading. To correctly determine the co-occurrence patterns, they suggest looking at the species pairs that have the highest CU to determine if they are segregating or aggregating (Ulrich and Gotelli, 2013). We looked at the top 2 percent of species pairs with the greatest absolute value of average CUs to determine their co-occurrence patterns (Appendix Fig. A1). Additionally, negative co-occurrence patterns can result from either species interactions or environmental heterogeneity, where species respond in dissimilar ways to underlying environmental variables. To determine the role of environment in structuring the negative species co-occurrences, we looked at the species pairs with the greatest absolute value of average CUs. The relative abundance of each species was plotted against the environmental gradient that they segregated along in the CCA (Appendix Fig. A1).
A total of 79 taxa were identified from the lakes and ponds in Wapusk National Park (Appendix Tables A1–A3). The most frequently found taxa were Conochilus spp. Daphnia tenebrosa, Heterocope septentrionalis, and Leptodiaptomus minutus. The frequency of occurrence for all taxa ranged from 1% to 77%, with an average of 22%. Richness was highest in rotifers and lower in cladoceran and copepod taxa (Fig. 2). Shannon-Weiner diversity was low, ranging from 0.06 to 1.02 with a mean diversity of 0.5 (Fig. 2). Presence/absence data for each lake/pond can be found in Appendix Table A4.
CCA allowed examination of the relationship between lake/pond environmental measures and variation in zooplankton community composition. The first and second axes represent 9.4% and 5.3% of the variation in zooplankton composition, respectively, and were both significant at p < 0.05. The significant environmental variables were lake area and chlorophyll a, which were mainly associated with axis 1, and pH, ionic composition, and TP, which were associated with axis 2 (Fig. 3). Lakes/ponds located with the three different habitats were separated along axis 1 (ANOVA, F = 17.5, p < 0.001), with lakes/ponds in the spruce forest region having low axis 1 scores, lakes/ponds in the coastal fen having high axis 1 scores, and lakes/ponds in the interior peatland being evenly distributed across the axis (Fig. 3). The spruce forest lakes/ponds have the lowest axis 1 scores (Tukey honest significant difference [HSD]: spruce forest — coastal fen, p < 0.001; spruce forest — interior peatland, p = 0.003) associated with large area, high chlorophyll a values, and species composition dominated by Keratella spp. Diaptomus nudus and Eucyclops serrulatus (Figs. 3 and 4). Lakes/ponds in the interior peatland region had intermediate axis 1 scores (Tukey HSD: interior peatland — coastal fen, p = 0.04; Fig. 3). Finally, the coastal fen lakes/ponds had the highest axis 1 scores associated with small area and low chlorophyll a values, dominated by the rotifers Trichocerca spp., Kellicottia spp., Synchaeta spp., and Notholca spp. and by harpacticoid copepods (Fig. 4).
Median and range of physical and environmental variables from the 92 lakes/ponds sampled.
Variation partitioning analysis revealed that spatial variables consistently explained more independent variation than environmental variables (Fig. 5 and Table 2). The environmental variables in this study explain less variation in zooplankton community composition than other studies of zooplankton distribution (Fig. 5). When investigating smaller spatial scales (coastal fen, interior peatland, and spruce forest) the independent effect of environment was rarely significantly different from 0, yet the independent effect of spatial variables was often significant (Table 2).
The co-occurrence analysis showed that there were significant negative co-occurrence patterns, or segregations, between zooplankton taxa. The standardized effect size of the co-occurrence analysis was highest in the rotifer assemblage (19.7) and lower in the cladoceran (6.3) and copepod (4.6) assemblages (Fig. 6). Most of the species pairs that had the greatest absolute value of checkerboard units (i.e., contributed the most to the standardized effect size) showed patterns of segregation and did not appear to be responding differentially to environmental gradients (Appendix Fig. A1).
The zooplankton communities of Wapusk National Park show comparable crustacean zooplankton diversity to other Subarctic/Arctic ponds at similar latitudes (Hebert and Hann, 1986). The number of crustacean zooplankton taxa we identified (42) represents an increase from the 25 taxa Hebert and Hann (1986) reported in Churchill. This increase is likely due to the larger spatial scale and increased diversity of habitats sampled. The water chemistry of the lakes/ponds (pH, conductivity, TP, TN, DOC) is within the range found in other shallow ponds in the Arctic, whereas the chlorophyll a concentrations were greater than most arctic sites (Rautio et al., 2011). The environmental gradients in the park are large, representing a shift from coastal habitat to inland boreal forest. While the data we collected and gradients we sampled were similar to or greater than those in other systems, where environment was able to explain large amounts of variation in zooplankton community composition (e.g., Rautio, 1998; Cottenie et al., 2001; Steiner, 2004; Strecker et al., 2008; Tavernini et al., 2009), our environmental variables explained relatively little of the variation in community composition.
Our initial goals involved determining the relative importance of environmental and spatial variables in structuring the zooplankton communities of Wapusk National Park. The three models that seemed most likely for our biota included dispersal limitation, species sorting, and mass effects. Of course, while it is possible that different species groups may not fit the same model (even in the same habitat) due to differing dispersal potential, we separated the three taxa in our analyses and found broadly similar results for all three. Species sorting entails a strong correlation between species composition and environmental conditions, and we found a small independent effect of environmental variables, especially compared to other studies of zooplankton distribution (Fig. 5). The spatially structured environment component of variation (ES) was relatively large compared to other studies (Fig. 5) and is likely caused by environmental differences between the three spatially separate habitat types (Fig. 3); however, if we consider the most extreme possibility that the entire ES component is an environmental signal, the environmental variables still explain less variation than almost all similar published studies (Fig. 5). There is potential for unmeasured environmental variables that vary with habitat type to increase the importance of the spatial factors (Cottenie et al., 2003), though given the extensive number of environmental factors we sampled, the importance of unmeasured variables is likely low. For example, although there were no data on predators in our analyses, only two of the lakes/ponds have known fish populations, therefore not having this data likely had minor impacts on the results. The independent spatial variables explained a large proportion of variation, suggesting that species sorting was not the sole mechanism occurring at the spatial scale we investigated. When considering the three regions separately, there is some evidence that the taxa are responding to environmental gradients due to high explanatory power of ES, but there is often a significant independent effect of space (Table 2).
The results of variation partitioning analysis; y is the exponent selected for use in the spatial model equation, # Variables represents the number of spatial variables selected for inclusion in the final model, E|S represents the independent variation explained by environmental variables, ES represents the shared variation explained by both environmental and spatial variables, S|E represents the proportion of independent variation explained by spatial variables, and U represents the proportion of unexplained variation.
Three metacommunity models that would involve a strong spatial signal are the “mass effects” metacommunity, where species are over-dispersed into habitats regardless of niche preferences; “neutral” metacommunities, where species are responding to dispersal constraints and stochastic competitive dynamics among equivalent species; and “dispersal limited” metacommunities, where species are unable to access ideal habitats due to dispersal constraints. Unfortunately, despite diametrically different causal mechanisms, these metacommunity models are notoriously difficult to distinguish in terms of pattern, and we have some evidence that argues in favor of each model. In support of the dispersal limitation paradigm, the amount of spatial variation explained decreased with the increasing ability for taxa to disperse (i.e., rotifers, with strong dispersal abilities [Caceres and Soluk, 2002], had spatial variation explaining less variance in composition than in the copepods, which are relatively weak dispersers [Frisch et al., 2012]). Furthermore, the negative co-occurrences among pairs of our species also support the hypothesis that the importance of space is due to dispersal limitation rather than true mass effects, as a mass effects scenario would suggest that species co-occur randomly. Conversely, other evidence suggests that mass effects are more likely than dispersal limitation. Most species had widespread distributions in the park, indicating that species likely have the potential to access all ponds. In addition, we also followed 23 ponds over three years (S. E. Arnott, unpublished data) and found that crustacean zooplankton species turnover averaged 28% per year, which is high compared to temperate lakes (Arnott et al., 1999) and suggests that dispersal (either through space or time via the egg bank) is high in Wapusk. In support of both the mass effects and the neutral model, the strong independent spatial signal could be evidence that environmental gradients are not providing strong selection against arriving species, although we acknowledge that this aspect of the models is difficult to assess in Wapusk because many of the environmental gradients were spatially structured. Perhaps the strongest argument in favor of a high dispersal interpretation of our data is the substantial evidence that zooplankton are capable of dispersing over large distances (De Meester et al., 2002), and the characteristics of our tundra-pond ecosystem (e.g., no trees, flat land, much wind, waterfowl, and ephemeral habitat) suggest that dispersal between lakes/ponds is high. Ultimately, however, we are unable to determine if the spatial signal we detected was due to low or high dispersal in this system; although, as we outline above, there seems to be more evidence for high dispersal. Recently, it has been argued that the models of mass effects, dispersal limitation, species sorting, and neutral dynamics fit on a continuum of dispersal and environmental heterogeneity, and it is difficult to categorize natural systems, as they often share characteristics with more than one paradigm (Logue et al., 2011; Winegardner et al., 2012). Furthermore, it is important to consider the interactions between the amount of dispersal, environmental gradients, and species interactions themselves.
Local interactions between species are also capable of having a strong effect on community composition (Diamond, 1975). Biotic resistance and/or priority effects have been extensively investigated experimentally (e.g., Shurin, 2000; Forrest and Arnott, 2006; Strecker and Arnott, 2010) and have a welldeveloped theoretical background (Steiner and Leibold, 2004). At the landscape level one way to determine if priority effects or local species interactions are important to community structure is to determine species co-occurrence patterns (Conner and Simberloff, 1979). Negative co-occurrence patterns are expected when species interactions or priority effects are excluding additional species from colonizing a habitat patch. As we found strong negative co-occurrence patterns for all three species groups across the landscape there is reason to expect that priority effects may be important in Wapusk. Conversely, a meta-analysis by Gotelli and McCabe (2002) showed that most biological communities show patterns of negative co-occurrence, however, these patterns may not necessarily be driven by competition, but rather by species responding differently to environmental gradients across the sampling scale. In our study, the species pairs that had the highest negative co-occurrence values (checkerboard units, CUs) were not separated along important environmental gradients, suggesting that different environmental tolerances were not driving the negative co-occurrence patterns and that other ecological processes (e.g., priority effects) were more important in causing the negative co-occurrences, as would be predicted by the neutral model. Indeed, neighboring lakes/ponds with similar chemistry in Wapusk have negatively co-occurring species pairs, indicating that we may even be seeing alternate stable equilibria. Finally, priority effects are expected to be particularly important in smaller habitats (e.g., ponds), as the populations are able to rapidly monopolize the location (Steiner and Leibold, 2004), and when the biota involved reproduce quickly and have the ability to rapidly colonize a space, for example, via a standing egg bank (De Meester et al., 2002). Altogether, these factors suggest that biotic interactions may play an important role in structuring communities, though experimental methods and temporal data will likely be needed to evaluate the importance of priority effects and dispersal in this system.
Overall, the three different habitat types had a large effect on zooplankton community composition, yet there were strong negative species co-occurrences suggesting that species interactions and priority effects may be important in determining species composition in this pond ecosystem. Given that we believe the system to have high levels of dispersal, we expect that zooplankton should be able to track changing environmental conditions. If the cover of coastal fen, interior peatland, and spruce forest habitats is affected by climate change, then this will likely have large implications for the distribution and abundance of aquatic taxa in Wapusk National Park.
We thank Justin Shead, Sheldon Kowalchuck, Heather Stewart, Brendan McEwan, Jill Larkin, Kevin Burke, David Walker, and Larry Gogal for assistance in the field. Thanks to Parks Canada, Manitoba Conservation, and Gogal Air for logistical support. Financial support was provided by an International Polar Year grant to Jon Sweetman, NSERC Discovery grant to Shelley Arnott, and a Polar Continental Shelf Project grant.
Physical data from each lake. All lakes are located in UTM zone 15N. Habitats are coded as CF for coastal fen, IP for interior peatland, and SF for spruce forest.
Chemical data from each lake, including: pH, conductivity (Cond, µS cm-1), total nitrogen (TN, mg L-1), total phosphorus (TP, mg L-1), temperature (Temp, C), dissolved oxygen (DO, mg L-1), alkalinity (Alk, mg L-1), chloride (Cl, mg L-1), sulfate (SO4, mg L-1), dissolved organic carbon (DOC, mg L-1), dissolved inorganic carbon (DIC, mg L-1), calcium (Ca, mg L-1), magnesium (Mg, mg L-1), potassium (K, mg L-1), sodium (Na, mg L-1), and silica (SiO2, mg L-1).
Chlorophyll-a data from each lake.