Ordination obscures the influence of environment on plankton metacommunity structure

The composition of plankton communities in individual habitats is often influenced by environmental conditions like pH or hydroperiod. At larger scales, environmental gradients can influence community structure across interconnected local communities. Detecting the role of environmental and spatial factors on metacommunity structure depends on the ordering of sites and species prior to analysis. We investigated this ordination in two wetland metacommunities; a well‐sampled, hyper‐diverse zooplankton metacommunity, and a Central American phytoplankton metacommunity. We calculated coherence, turnover, and boundary clumping to classify the structure of the metacommunity, and we propose a statistic that responds to variation in both coherence and turnover. Traditional ordination approaches failed to discern metacommunity structure, while significant structure existed along abiotic gradients in both zooplankton and phytoplankton systems. This shows that abiotic controls on community composition may not be detectable with traditional analyses, and suggests an alternative of ordering sites by known abiotic gradients.


Scientific Significance Statement
Freshwater systems are often thought to operate as metacommunities, or discrete communities connected through dispersal. There are different ways to analyze metacommunities to understand the underlying factors controlling structure, most of which use some form of ordination. Here, we modify one analysis, showing where the effect of environmental controls on metacommunity structure in plankton metacommunties can be obscured by traditional ordination and we propose a novel way to better visualize metacommunity structure. of metacommunity ecology. This combination of theory from community ecology and biogeography has been applied to a range of complex ecosystems, including zooplankton communities in ponds and lakes (Cottenie et al. 2003), parasite communities of host species (Dallas and Presley 2014), and microbial communities of pitcher plants (Miller and Kneitel 2005). While conceptually compelling, there is not agreement about what quantitative framework best captures the factors influencing metacommunity structure (i.e., quantifiable patterns of nestedness or clustering of patches in environmental or geographic space) (Ulrich and Gotelli 2013).
One approach to quantifying metacommunity structure is the pattern-based framework proposed by Leibold and Mikkelson (2002). This framework characterizes metacommunity structure according to three statistics: coherence, turnover, and boundary clumping. Coherence measures the continuity of a species range, quantified as the number of absences within the species range boundaries (i.e., embedded absences). Significantly positive coherence (i.e., low number of embedded absences) is a necessary condition for the other two statistics to be interpretable in the Leibold and Mikkelson (2002) framework. The second statistic, turnover, measures the tendency of species to replace one another at their range boundaries, and is closely related to measures of species co-occurrence [i.e., checkerboard score; Stone and Roberts (1990)]. The final statistic of the Leibold and Mikkelson (2002) framework is boundary clumping, which measures the tendency of species range boundaries to occur together. Inferences may be made by comparing each of these statistics (apart from boundary clumping) to a null distribution generated by permutation of species occurrences subject to metacommunity-level constraints. The metacommunity statistics (coherence, turnover, and boundary clumping) are calculated on the empirical interaction matrix, compared to the null distributions of statistics from permutation analysis, and the metacommunity is classified as one of six (Leibold andMikkelson 2002) [or ten, Presley et al. (2010)] possible structures. This framework has been criticized, however, due to variation in performance on simulated data and the known sensitivity of estimated statistics to matrix ordering (Ulrich and Gotelli 2013). This critique could be addressed by using statistics that are invariant to matrix order [e.g., Barber's modularity Q; Barber (2007)], or by ordering sites and species along known environmental gradients (Gotelli and Ulrich 2012). Despite these potential shortcomings, we believe the Leibold and Mikkelson (2002) framework may nevertheless be of considerable value. Particularly, it has the virtue of providing quantities that map clearly to basic metacommunity concepts and may be modified and extended as needed to represent new ideas (Presley et al. 2010).
We used the Leibold and Mikkelson (2002) framework, and addressed the critique of Ulrich and Gotelli (2013) in two novel ways. First, we compared metacommunity classifications based on the reciprocal averaging approach of Leibold and Mikkelson (2002) to an approach that orders sites and species along known environmental gradients. Reciprocal averaging (also known in community ecology as "correspondence analysis") places sites with similar species, and species with similar distributions, closer together in the site-by-species interaction matrix. Commonly, these ordination scores are related to environmental or spatial covariates to infer the relative influence of covariates on metacommunity structure (Willig et al. 2011;L opez-Gonz alez et al. 2012;Dallas and Presley 2014). However, correlations between covariates and ordination scores may be misleading, while ordering sites by known environmental gradients may lead to different conclusions. Specifically, assessing metacommunity structure along multiple environmental gradients could result in different classifications (e.g., random, nested), which could suggest conflicting influences of environmental covariates on metacommunity structure. We demonstrate the sensitivity of the analysis to matrix ordering, and recommend ways to address this issue. Second, we extend the Leibold and Mikkelson (2002) framework by advocating the use of a continuous measures of the magnitude of coherence and turnover to quantify variation in metacommunity classifications within a continuous space [see Heino et al. (2015)]. In our analysis, a metacommunity occupies a single point within a continuous space of coherence and turnover. Successional changes, species extinctions, immigration events, and species invasions may not alter metacommunity classification, but will affect the location of the metacommunity in the coherenceturnover phase space. Previous studies have visualized different metacommunities in this phase space (Heino et al. 2015).
We propose instead to study the location in this space of a single metacommunity, ordered according different environmental variables. The distance between points in coherenceturnover space offers a novel way to examine the influence of time or null model choice on metacommunity structure. We examined two aquatic metacommunity systems (one zooplankton and one phytoplankton) to investigate our extension to the Leibold and Mikkelson (2002) framework. Plankton metacommunities are well-studied, although the relative importance of spatial (Shurin 2000;Rojo et al. 2016), and environmental (Cottenie et al. 2001(Cottenie et al. , 2003 covariates remains unclear (Soininen et al. 2007;Dallas and Drake 2014). Both phytoplankton and zooplankton community data were obtained through repeated sampling of two series of wetlands. Zooplankton were sampled from 14 Carolina Bays, and phytoplankton from 30 wetlands in Costa Rica and Nicaragua (Rojo et al. 2016). For both data sources, we analyzed metacommunities ordered by traditional ordination analysis, and compared the outcome when sites were ordered along known environmental gradients. We demonstrated that (1) classification of plankton metacommunity structure is sensitive to the order of sites and species, and (2) plankton assemblages were structured by environmental covariates that were only weakly related to ordination scores. Taken together, these findings provide evidence of environmental structure in zooplankton and phytoplankton communities, reveal the pitfalls of the traditional ordination approach, and demonstrate a useful extension to the Leibold and Mikkelson (2002) framework for characterizing metacommunity structure.

Zooplankton communities
Zooplankton communities were sampled every other week between January 2009 and February 2011 from 14 ephemeral water bodies within the Savannah River Site (SRS), a nuclear reserve owned by the United States Department of Energy located in the upper coastal plain of South Carolina (Fig. 1). These ephemeral wetlands, also referred to as Carolina Bays, are an ideal system to examine metacommunity structure, as communities are extremely diverse, well-sampled, and occupy discrete habitats. These ephemeral wetlands usually have no natural surface drainage and depend largely on precipitation for filling, resulting in typically low conductivity and acidic conditions (Newman and Schalles 1990). Sampled wetlands (n 5 14) are typically inundated in the winter and dry in summer [although some stay inundated year round except in drought years; Sharitz (2003)]. These wetlands are strikingly species-rich, as single wetlands can have as many as 60 zooplankton species (Zokan 2016), and over 100 zooplankton species recorded throughout the system of ephemeral wetlands (DeBiase and Taylor 2005). A total of 84 zooplankton species were recorded over the course of 49 potential sampling events over the time period examined here. Further information related to water quality and sampling protocols can be found in the Supporting Information.

Phytoplankton communities
Phytoplankton communities were sampled by Rojo et al. (2016). Sites were sampled at three times between September 2010 and June 2011, corresponding to periods of inundation and dessication of wetlands. Here, we combine data from two sampling periods (January and June 2011), excluding the September 2010 sampling as some of the wetlands sampled later were not sampled in September. This resulted in a total of 295 phytoplankton species across the 30 sampled wetlands. See Rojo et al. (2016) for more detailed information on sampled sites and methodology, as well as experimental data [provided as Supporting Information in Rojo et al. (2016)].

Metacommunity analyses
We used the Leibold and Mikkelson (2002) framework to analyze zooplankton species assemblages among ephemeral bays, and phytoplankton assemblages among wetlands, both represented as site-species matrices of wetlands (rows) and species (columns). The ordering of sites and species can influence statistical values of the Leibold and Mikkelson (2002) framework, and may alter metacommunity classification. Sites were ordered by either reciprocal averaging scores or observed environmental gradients. Reciprocal averaging concentrates species occurrences along the matrix diagonal to group species with similar ranges together, and sites with similar species assemblages together. Column ordination was based on reciprocal averaging in all cases, while row ordering was based on central tendencies of measured environmental variables. For zooplankton communities, environmental variables included mean values of conductivity, water depth, hydroperiod (length of time the bay contains water), pH, and temperature. For phytoplankton communities, sites were ordered based on mean values of conductivity, water depth, pH, temperature, dissolved oxygen, chlorophyll a, bicarbonate, nitrate, and phosphate. Following Leibold and Mikkelson (2002), three statistics were calculated on each ordered site-species matrix; coherence, turnover, and boundary clumping. Coherence is quantified as the number of absences between the extremes of species ranges, and measures the tendency for species ranges to be constrained to a subset of similar sites (Leibold and Mikkelson 2002). Turnover is a measure of overlap at the extremes of species ranges, and is calculated after species ranges are made completely coherent (embedded absences are removed). Boundary clumping measures the tendency of the extremes of species ranges to coincide with the range extremes of other species, and is quantified using Morisita's index (Leibold and Mikkelson 2002). This measure is conceptually similar to the field of community detection and measures of modularity in graph theory, although measures of modularity are typically invariant to site-species matrix ordering (Barber 2007). Here, we use Morisita's index since we aimed to examine the sensitivity of the Leibold and Mikkelson (2002) framework to matrix ordering.
Statistical significance was determined relative to the distribution of the statistic calculated on randomized matrices (i.e., a null distribution) for coherence and turnover. Significance of boundary clumping (Morisita's index) was assessed relative to a chi-squared distribution. The randomization of species occurrences, much like the ordering of the interaction matrix, can influence statistical results and metacommunity classification. We chose a conservative null model [sequential swap algorithm; Gotelli and Entsminger (2003)] that maintains the number of sites occupied for each species and the number of species in each local community (also called a fixed-fixed null model). However, we relax this null model in the Supporting Information, and instead assign species occurrences proportional to the number of observed occurrences [called a fixed-proportional null model; Gotelli (2000)].
Significance was determined by calculating each statistic for a set of 1000 null matrices, and comparing the statistic calculated on the empirical matrix using a z-test. However, significance does not address the magnitude of divergence from the null expectation. To do so, we used the z statistic as a measure of the divergence of our empirical statistic from the statistic's null distribution. We computed the z statistic as the mean simulated statistic minus the observed statistic divided by the standard deviation of the simulated statistic.

Relationship of environmental gradients to reciprocal averaging gradient
Values for each environmental gradient were compared to ordination scores (first axis of reciprocal averaging ordination) to determine the strength of association between Null values of c (denoted with an upper-case C) were calculated based on the swap null model, which fixes row and column totals. These null values (mean 5 C , variance 5 r C ) were compared to the empirical c to assess significance and calculate divergence from null expectation (z). "Recip. avg." refers to the metacommunity ordered based on reciprocal averaging ordination scores. Null values for b (denoted with a B) were calculated based on the swap null model, which fixes row and column totals. These null values (mean-5 B, variance 5 r B ) were compared to the empirical b to assess significance and calculate divergence from null expectation (z). "Recip. avg." refers to the metacommunity ordered based on reciprocal averaging ordination scores. environmental gradient and the latent ordination gradient. The underlying assumption made by previous studies (Willig et al. 2011;L opez-Gonz alez et al. 2012) is that environmental factors important for metacommunity structure would be strongly correlated with the ordination scores. To test this, we examined correlations between environmental gradients and the first ordination axis, hypothesizing that environmental gradients along which the metacommunity was strongly structured would be more strongly associated with the structuring gradient obtained from the reciprocal averaging ordination.

Metacommunity structure along environmental gradients
When the site-species matrix was ordered by reciprocal averaging ordination scores, the zooplankton metacommunity did not differ from the expectation under the fixed-fixed (sequential swap) null model with respect to the number of embedded absences in species ranges (i.e., coherence; Table 1) or species turnover (Table 2). Based on the framework of Leibold and Mikkelson (2002), this metacommunity would be considered random. However, when sites were ordered based on pH, the metacommunity contained fewer embedded absences than expected by chance, indicating non-random metacommunity structure (Table 1).
Apart from being non-random, the metacommunity ordered based on pH had significantly positive turnover (more species turnover than expected under the fixed row and column null model; Table 2) and boundary clumping (pH M 5 10.98, df 5 84, p < 0.0001) suggesting that species replaced one another in discrete groups across the pH gradient (i.e., a Clementsian metacommunity). Results for the fixed-proportional null model were similar, although the metacommunity ordered along the reciprocal averaging gradient and conductivity become non-random ( Fig. 3 and Supporting Information).
Similar results were found for the phytoplankton metacommunity, where the number of embedded absences (i.e., coherence) in the traditionally ordered metacommunity was not significantly different than expected under the fixed-fixed null model, while the same metacommunity ordered along gradients of dissolved oxygen and pH yielded significantly non-random metacommunities (Table 1). Moreso, metacommunities ordered based on dissolved oxygen and pH had significantly greater turnover than expected (Table 2) and significantly clumped range boundaries (dissolved oxygen M 5 9.98, p < 0.001; pH M 5 12.09, p < 0.001). When sites and species were ordered using the traditional reciprocal averaging ordination, the metacommunity did not differ from the null expectation with respect to coherence, and would be classified as random. Thus, the same metacommunity ordered along the traditional reciprocal averaging gradient and two environmental gradients (pH and dissolved oxygen), produce two starkly different results.

Magnitude measures
The two metacommunity statistics that were compared relative to null models (i.e., coherence and turnover) allow for the creation of quantitative measures for the magnitude of effect, as described above. These measures create a phase space, where metacommunities ordered by different factors occupy different regions (Fig. 2). This continuous space adds to the simple classification of metacommunities by providing information on the distance between classification states. For instance, the zooplankton metacommunity was Clementsian when ordered by pH, suggesting that species replaced one another in discrete communities across the pH gradient, Fig. 2. Zooplankton (a) and phytoplankton (b) metacommunities were structured differently when sites were ordered by environmental variables or reciprocal averaging gradient (red point; "Recip. avg."). The magnitude of the difference between the null expectation and the empirical statistics for coherence (y-axis), and turnover (x-axis) determines the location of the metacommunity in the phase space. Gray lines separate regions based on statistical significance (a < 0.05) relative to a null expectation, where larger positive values of coherence magnitude indicate fewer embedded absences relative to the null model, and larger negative turnover magnitude indicates more species replacements than expected under the null model. Fig. 3. The choice of null model influenced the results of both zooplankton (a) and phytoplankton (b) metacommunities. Results from both null models are plotted in the coherence and turnover phase space, with dark gray lines connecting the same environmental gradients for both fixed-fixed (gray triangles) and fixed-proportional (blue circles) null models. Gray horizontal and vertical lines correspond to significance levels at a 5 0.05.
which can be visualized in phase space (Fig. 2). The position of a metacommunity in this phase space changes when using a different null model (Fig. 3), providing a means to evaluate the robustness of metacommunity classifications.

Associations between environmental covariates and the reciprocal averaging gradient
The environmental gradient responsible for structuring the zooplankton metacommunity was pH, as discussed above (Tables 1, 2). However, the metacommunity was classified as random when ordered using the traditional analytical approach (i.e., ordering by reciprocal averaging scores). The gradient obtained from reciprocal averaging was unrelated to the conductivity gradient ( Fig. 4; t 5 0.96, df 5 12, p 5 0.356), and only weakly correlated with the pH gradient ( Fig. 4; t 5 22.23, df 5 12, p 5 0.046). Hydroperiod did not influence zooplankton assemblages, but was strongly related to other environmental covariates (i.e., mean bay depth and temperature).
Similarly for the phytoplankton metacommunity, the reciprocal averaging gradient did not result in the same metacommunity classification compared to ordering sites along known environmental gradients. The phytoplankton metacommunity would be considered random when ordered along the traditional ordination gradient. Meanwhile, environmental gradients of dissolved oxygen and pH resulted in significantly structured metacommunites. Both dissolved oxygen (rho 5 0.07, df 5 28, p 5 0.03) and pH (rho 5 0.40, df 5 28, p 5 0.03) were weakly related to the reciprocal averaging ordination gradient (Fig. 4).

Discussion
The zooplankton metacommunity was classified as Clementsian along a gradient of pH, suggesting that discrete subsets of the zooplankton community have overlapping pH tolerances, resulting in a modular interaction pattern. However, the metacommunity was classified as random when sites were ordered by their reciprocal averaging scores [i.e., the traditional Leibold and Mikkelson (2002) approach]. Further, the phytoplankton metacommunity was classified as Clementsian when sites were ordered along gradients of dissolved oxygen and pH, while the metacommunity was classified as random when sites were ordered along the reciprocal averaging gradient. Regardless of metacommunity classification, the structuring environmental gradients for both zooplankton and phytoplankton were weakly related or unrelated to the reciprocal averaging gradient, and correlations often differed in the direction of the relationship. This suggests that the reciprocal averaging gradient was able to capture some environmental variation, but still masked the influence of individual environmental gradients. Perhaps more importantly, this suggests that the reciprocal averaging gradient may not capture all of the relevant biological information, leading to incorrect conclusions about metacommunity structure.
The Leibold and Mikkelson (2002) framework has been criticized as too simplistic, and potentially inaccurate, based on a simulation study examining several commonly used measures of species aggregation and segregation (Ulrich and Gotelli 2013). However, Ulrich and Gotelli (2013) found that the Leibold and Mikkelson (2002) statistics performed comparably to other measures of aggregation and segregation, despite being taken out of the context of the analytical chain proposed in Leibold and Mikkelson (2002). Specifically, matrices generated randomly, including those that would be classified as random, were used to assess the statistical properties of the other statistics in the Leibold and Mikkelson (2002) framework. Here, we alleviate some of the concerns raised in Ulrich and Gotelli (2013) by extending the Leibold and Mikkelson (2002) framework; introducing a continuous measure of coherence and turnover, calculating statistics relative to a number of realistic null models (see Supporting Information), and ordering sites along known biological gradients (Gotelli and Ulrich 2012).
Using our Leibold and Mikkelson (2002) extension, we found pH is important to the structure of ephemeral pond zooplankton communities. Previous work in lakes in the northeast United States (Dallas and Presley 2014), boreal shield lakes (Derry et al. 2009), and interconnected ponds (Cottenie et al. 2001) has suggested the importance of pH as a determinant of zooplankton community structure. This is arguably because zooplankton species have relatively narrow regions of pH conditions at which they can exist (Holt et al. 2003). The Carolina bays studied are rain-driven, and are typically more acidic than more permanent water bodies that may have a groundwater source. Hydroperiod can influence zooplankton diversity (Serrano and Fahd 2005;Zokan and Drake 2015), resilience (Angeler and Moreno 2007), and colonization dynamics (Chaparro et al. 2016) in ephemeral waterbodies. We failed to find an effect of hydroperiod on zooplankton community structure in our system. This may be because we didn't explicitly examine temporal patterns in zooplankton community composition (i.e., succession). To date, few studies of metacommunities have considered how metacommunity structure may change over time (but see Keith et al. 2011;Newton et al. 2012). However, the Leibold and Mikkelson (2002) framework may be extended to explain successional metacommunities by examining metacommunities as multi-layer networks (Pilosof et al. 2015), or examining changes in the coherence-turnover phase space we suggest here as a measure of metacommunity change through time.