Community stability:
We examined the temporal dynamics across six years (2010-2015) of the six coastal pond fish and invertebrate communities in ordination space using non-metric multidimensional scaling (nMDS). Testing for statistical differences in community composition across pond and year, we used a permutational multivariate analysis of variance (PERMANOVA). We used the functions 'metaMDS' and 'adonis' of the vegan R package to produce and plot the nMDS and calculate PERMANOVA respectively (R Development Core Team 2013). To calculate community instability, we examined variation in ordination space, with 2D Euclidean distance representing the change in the community composition (identity and abundance) across years (Clarke 1993, Kroeker et al. 2013, Lamothe et al. 2019). We adapted the three metrics of community instability from Mellin et al. (2010): (1) convex hull area, (2) average year-to-year distance, and (3) average year-to-centroid distance. The first metric, convex hull area, is the smallest possible convex polygon drawn in ordination space that encompasses all communities through time for each pond. Convex hull area gives one overarching instability value for each pond and thus quantifies the total change in community composition over time. The second metric, year-to-year distance, is the average year-to-year Euclidean distance for each pond community. Because this metric tracks changes in community composition from year-to-year, small annual shifts in the trajectory of the community due to incremental but directional environmental change will yield a low instability value, even if these small annual shifts yield a large change in community composition across all years. Hence, this metric quantifies the degree of incremental change in community composition over time. The third instability metric is the average year-to-centroid distance, which is the average distance between annual community structure and the mean community structure across all years (centroid). This metric thus tracks how much a community deviates annually from its average composition across all years.

After calculating the instability metrics, we then examined the relationships between pond size or ocean distance and each metric. We fit linear regressions to the pond size-/ocean distance-convex hull instability relationship. Because both year-to-year and year-to-centroid distance produced estimates for each pond in each year, we fit an ANCOVA to determine the strength and the persistence of the relationship between these instability metrics and both pond size and ocean distance (covariates) across years (factor). To confirm that data are not correlated across sampling events, we ran a Durbin-Watson temporal autocorrelation test as well as plotted the residuals of the linear models using the 'acf' function (Appendix S3 of Yeager et al., 2020). For all instability analysis, both pond size and ocean distance were log transformed.

Path analysis:
To understand whether community instability was directly affected by pond size and ocean distance or whether it was mediated through the relationship with species richness, we conducted path analysis using the lavaan R package (R Development Core Team 2013). Path analysis is a procedure that estimates the direct and indirect relationships between the exogenous variables pond size and ocean distance and the endogenous variables species richness and instability by computing standardized (path) coefficients using a series of multiple linear regressions (Grace and Bollen 2005, Grace et al. 2012).

We conducted two path analyses based on two of the three calculated instability metrics: (1) year-to-year distance and (2) year-to-centroid distance. Due to their intrinsic properties, these metrics provided different information about community instability. Year-to-year distance describes how instability changes incrementally with time, whereas year-to-centroid distance represents total instability across all years. For each path analysis, we calculated the direct effects of both pond size and ocean distance on instability. The indirect effects of pond size and ocean distance on instability were obtained by computing the product of the standardized coefficient linking pond size or ocean distance to species richness and the standard coefficient linking species richness to instability. To understand the degree to which the relationship between pond properties and instability was driven by direct or indirect effects, we calculated the total effect of pond size and ocean distance by summing their direct and indirect effects. For both path analyses, we determined whether the observed and the predicted covariance matrices differed (Appendix S4 of Yeager et al., 2020) by computing the summary  statistic and its statistical significance based on 1,000 iterations using the Yuan bootstrap method (Yuan et al. 2007). Additionally, we ran the path analyses using the rarefied species richness values to ensure that our results were not a product of variation in the number of sampling stations across ponds (Appendix S2 of Yeager et al., 2020).