Introduction

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

Methods

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

Results

Lakes and ponds were situated across costal fen, interior peatland, and spruce forest (Fig. 1). Among the 92 lakes/ponds the geographical and environmental conditions varied considerably (Table 1).

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

TABLE 1

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

Discussion

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.

FIGURE 1.

Map of the study area, Wapusk National Park. Squares represent lakes/ponds in the spruce forest region, triangles represent lakes/ponds in the interior peatland region and circles represent lakes/ponds in the coastal fen region.

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

FIGURE 2.

The distribution of (a) total community Shannon-Weiner diversity, (b) total richness, (c) rotifer genera richness, (d) copepod species richness, and (e) cladoceran species richness in the 92 lakes/ponds.

FIGURE 3.

CCA biplot of site scores for the 92 ponds. λ represents the percentage of variation explained by each axis. Squares represent lakes/ponds in the spruce forest region, stars represent lakes/ponds in the interior peatland region, and open circles represent lakes/ponds in the coastal fen region.

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.

FIGURE 4.

CCA biplot of species scores for the 92 lakes/ponds. λ represents the percentage of variation explained by each axis. Only species with scores > 0.7 were plotted for clarity. Abbreviations: L.SET: Latona setifera; L.TYR: Leptodiaptomus tyrelli; S.CRY: Sida crystalline; BOS: Bosmina spp.; KEL: Kellicottia spp.; SYN: Synchaeta spp.; HAR: Harpacticoid copepods; TRI: Trichocerca spp.; NOT: Nothlca spp.; P.PED: Polyphemus pediculus; E.SER: Eucyclops serrulatus; KER: Keratella spp.; D.NUD: Diaptomus nudus.

FIGURE 5.

Results of variation partitioning analysis from studies of crustacean (i.e., copepod and cladoceran) zooplankton distribution. Black regions represent the independent effect of environment, the hatched regions represent the shared variation between environmental and spatial variables, the shaded regions represent the independent effect of spatial variables, and the white regions represent unexplained variation in community composition.

FIGURE 6.

Results of co-occurrence analysis for different zooplankton taxa. Significant negative co-occurrence patterns are detected at standardized effect size values above 1.96 (dashed line) at p < 0.05.

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.