Land Use and Salinity Drive Changes in SAV Abundance and Community Composition

Conserving and restoring submerged aquatic vegetation (SAV) are key management goals for estuaries worldwide because SAV integrates many aspects of water quality and provides a wide range of ecosystem services. Management strategies are typically focused on aggregated abundance of several SAV species, because species cannot be easily distinguished in remotely sensed data. Human land use and shoreline alteration have been shown to negatively impact SAV abundance, but the effects have varied with study, spatial scale, and location. The differences in reported effects may be partly due to the focus on abundance, which overlooks within-community and among-community dynamics that generate total SAV abundance. We analyzed long-term SAV aerial survey data (1984–2009) and ground observations of community composition (1984–2012) in subestuaries of Chesapeake Bay to integrate variations in abundance with differences in community composition. We identified five communities (mixed freshwater, milfoil-Zannichellia, mixed mesohaline, Zannichellia, and Ruppia-Zostera). Temporal variations in SAV abundance were more strongly related to community identity than to terrestrial stressors, and responses to stressors differed among communities and among species. In one fifth of the subestuaries, the community identity changed during the study, and the probability of such a change was positively related to the prevalence of riprapped shoreline in the subestuary. Mixed freshwater communities had the highest rates of recovery, and this may have been driven by Hydrilla verticillata, which was the single best predictor of SAV recovery rate. Additional species-specific and community-specific research will likely yield better understanding of the factors affecting community identity and SAV abundance, more accurate predictive models, and more effective management strategies.


Introduction
Submerged aquatic vegetation (SAV) is a foundational element in estuarine and shallow marine environments . The physical, biological, and biogeochemical impacts of SAV on aquatic environments shape ecosystem dynamics, including carbon sequestration (Fourqurean et al. 2012), water clarity (Gruber et al. 2011), primary and secondary production, and food web structure (Heck et al. 2003;Larkum et al. 2006). SAV is globally imperiled because its need for light and adequate substrate makes it sensitive to anthropogenic impacts . The sensitivity, ecological importance, and economic value of SAV together make its conservation and restoration a common target for monitoring and management programs in aquatic ecosystems worldwide Waycott et al. 2009).
Land use change and shoreline alteration can negatively impact SAV by changing the movement of sediment and nutrients from the land to the water (Bilkovic et al. 2006;Li et al. 2007;Blake et al. 2014). The effects of these stressors on SAV seem dependent on spatial scale and interaction with other Communicated by Masahiro Nakaoka environmental conditions (Patrick et al. 2014;Patrick et al. 2016). In Chesapeake Bay, shoreline armoring was observed to have stronger negative relationships with SAVabundance (a metric of SAV density and areal coverage; see BMethods^for more detail) in more saline waters and was more apparent in watersheds with less human land use. Negative relationships between the proportions of developed and agricultural land and overall SAV abundance were statistically significant but had weak explanatory power, suggesting that other factors also affect the linkages between human activities and SAV (Patrick et al. 2014;Patrick et al. 2016).
Incorporating temporal variability into analyses may help yield better understanding of stressor-response relationships than reported for analyses of temporally averaged SAV abundance (Patrick et al. 2014;Patrick et al. 2016). Historical impacts can persist in systems as legacy effects, obscuring or masking current relationships (Foster et al. 2003;Allan 2004). SAVs are hypothesized to exhibit positive density dependence (van der Heide et al. 2011), because larger beds can buffer against changes in water quality (Gruber et al. 2011;Gurbisz and Kemp 2014;Gurbisz et al. 2016), so areas with large SAV beds may be buffered against increases in development and able to maintain SAV populations. In contrast, areas that have lost much of their SAV to an episodic disturbance such as a hurricane may be unable to recover if chronic anthropogenic stressors are present on the surrounding landscape or if these locations are recruitment limited (Kendrick et al. 2012).
Similarly, incorporating species and community information into SAV models may clarify confusing relationships reported between terrestrial stressors and total SAVabundance (Kemp et al. 2004;Patrick and Weller 2015). Land use is hypothesized to impact SAV when it promotes suspended sediments and eutrophication that reduce light availability or when shoreline land use reduces potential habitat by changing nearshore wave energy and eroding shallow water habitat (Koch 2002;Strayer and Findlay 2010;Gittman et al. 2015). Species-specific differences in habitat needs may affect the strength of relationships between anthropogenic stressors and SAV abundance (Batiuk 2000;Kemp et al. 2004;Patrick and Weller 2015). Stressors can also eliminate sensitive species so that species composition shifts to tolerant or invasive species even if coverage does not change. We need to understand these species-specific responses to understand and effectively manage the effects of human activities on SAV.
Chesapeake Bay is an excellent study system to explore how community composition influences the relationships of terrestrial anthropogenic stressors with spatial and temporal variations in SAV abundance. Bay-wide SAV abundance has been mapped annually since 1984, and prior research informs our understanding the system (Orth et al. 2010;Ruhl and Rybicki 2010;Williams et al. 2010;Gurbisz and Kemp 2014;Patrick et al. 2014). The many tributary subestuaries to the bay provide replicate study units with strong differences in human land use (agricultural or urban land cover or shoreline armoring), and different regions have different temporal trends and interannual fluctuations (Orth et al. 2010;Patrick and Weller 2015). The rich spatiotemporal data and the wide range of stressor conditions facilitate testing how stressors interact with communities through time. Moore et al. (2000) identified four major SAV communities in Chesapeake Bay: a Zostera marina community, a Ruppia maritima community, a Potamogeton spp. community, and a mixed freshwater community. They were comprised of 23 different species that differ in growth morphology, substrate needs, sensitivity to changes in light availability, and other characteristics (Batiuk 2000). For example, some species such as Hydrilla verticillata are able to form dense canopies at the water's surface (Patrick and Weller 2015). Models relating human stressors to SAV have been improved by allowing stressor-response relationships to differ among salinity zones (Li et al. 2007;Patrick et al. 2016), and considering communities rather than salinity zones may yield further improvement (Patrick and Weller 2015).
We incorporated temporal dynamics and SAV communities into one series of analyses designed to determine how terrestrial human stressors affect SAV abundance in estuaries. The analyses examined Chesapeake Bay SAV between 1984 and 2009, a period when SAV showed dynamic changes following the catastrophic collapse in the 1970s. We focused on determining how land cover, land cover change, and shoreline armoring related to changes in SAV abundance within subestuaries and how these relationships are influenced by the species and communities present within the subestuaries. We also investigated how spatial patterns in stressors related to temporal patterns of turnover in SAV community identity.
First, we hypothesized that the strength of terrestrial anthropogenic stressors (such as the proportion of agricultural land use) in the local watershed is negatively related to the SAV abundance in the receiving subestuary. Subestuaries with stronger terrestrial stress (such as more agricultural land) will have slower rates of SAV recovery than those with less stress, and subestuaries where terrestrial stressors are increasing will have slower recovery than subestuaries where stressors are not increasing. Second, we hypothesized that subestuaries in which the SAV community is dominated by species with lower light requirements, broader substrate tolerances, or canopy forming growth morphologies will be more resilient to stressors than subestuaries with the opposite characteristics. Finally, we hypothesized that subestuaries with stronger terrestrial stressors will have lower community composition fidelity, because disturbance opens habitat for the invasion and proliferation of stressor tolerant exotic SAV species such as H. verticillata and Myriophyllum spicatum.

Study Region and Subestuaries
The Chesapeake Bay has 18,804 km of shoreline and a 164,000 km 2 drainage basin. Salinity ranges from tidal freshwater to 25 near the Atlantic Ocean outlet. Our study focused on 95 subestuaries-relatively small embayments (median size 11.8 km 2 ) connected to the main channel of the bay or to one of the major tributary rivers (Fig. 1). Each subestuary has its own local watershed, and among the subestuaries, local watershed land cover ranges widely between 0 and 91% developed land and 0.20-64% cropland (Patrick et al. 2014). Salinity zone was assigned to each subestuary ( Fig. 1) from the Chesapeake Bay Program (CBP) segmentation scheme (2004), and mean salinity was estimated from the CBP water quality monitoring program data (Patrick and Weller 2015). Additional subestuary characteristics were summarized from bathymetric maps (mean depth in m and volume in m 3 ), shoreline maps (shoreline length in km, subestuary mouth width in km, shoreline fractal dimension, subestuary area in km 2 , watershed/estuary area ratio), and long-term weather data (mean annual precipitation in mm). Mean tidal range and relative sea level were assigned to each subestuary from an interpolation of tidal range data from 657 tide stations and an interpolation of changes in annual mean water elevations over time measured at 28 National Ocean Service data stations. Details on the data sources and summary calculations are published elsewhere (Li et al. 2007;Patrick et al. 2014).

Terrestrial Stressors
We summarized published digital maps to characterize each subestuary and its local watershed. Land cover percentages (developed land, cropland, forested land, and wetland) were calculated for the local watershed of each subestuary in 1984, 1992, 2001, and 2006. Land cover was summarized from the Chesapeake Bay Watershed Land Cover Data Series which were derived from Landsat 5 Thematic Mapper and Landsat 7 Enhanced Thematic Mapper satellite imagery (Irani and Claggett 2010). The layers contain 16 land use and land cover classes which were further summarized as follows: Total development was the sum of open, low, medium, and high development classes. Total forest was the sum of deciduous, evergreen, woody wetland, and scrub shrub. Subestuaries were also classified by the dominant land cover regime ( Fig. 1) in their local watersheds into three categories: forested (>60% forest), human impacted (>50% developed or >40% cropland), or mixed (not fitting in the first two categories) (see Patrick et al. 2015). Shoreline alteration information for each subestuary was derived from the digital shoreline situation maps in the Comprehensive Coastal Inventory Program (http://www.vims.edu/ccrml). For each subestuary, we summarized the % shoreline which was hardened (bulkhead, riprap, or other armoring structures) and the density of built structures along the shoreline (e.g., docks or boat ramps in structures/km of shoreline).

SAV Abundance
The abundance of submerged aquatic vegetation coverage in each subestuary was summarized annually from 1984 to 2009 digital maps derived from aerial photography (http://www. vims.edu/bio/sav). The maps outline areas occupied by SAV and provide the density of SAV in each mapped bed (five density classes with coverage estimated as 0, 0-10, 10-40, 40-70, or 70-100). To quantify overall SAV abundance in subestuaries, we summarized the density weighted occupied habitat (hereafter referred to as abundance) for each subestuary for each year from 1984 to 2009 (except 1988 when no data were collected). Abundance was calculated as the sum of the area of each SAV bed multiplied by the midpoint of the density category assigned to the bed, then divided by the potential habitat area in the subestuary (see Li et al. 2007). This calculation normalizes for differences among subestuaries in size or in the area of potential habitat available. Potential habitat was estimated as the area of submerged habitat <2 m deep, in a bathymetric map (Cohen 1994), minus areas designated as no-grow zones by the CBP (USEPA 2003) plus any areas less than 2 m deep or within the no-grow zones where SAV has actually been observed (details in Patrick et al. 2014).

Community Classification
We summarized a long-term database of SAV species observations to characterize the SAV communities. The species identifications came from VIMS SAV ground survey data (http:// web.vims.edu/bio/sav/field_observations.html), which are collected annually by researchers, natural resource managers, and trained volunteers. The observers report species identifications coded by location and date, and the data from 1984 to 1995 have been previously used to characterize the communities in SAV beds across Chesapeake Bay (Moore et al. 2000). To characterize the average species composition within each subestuary, we combined all of the ground survey observations from 1984 to 2012 to create a summed subestuary × species occurrence observation matrix for the entire period of record. The number of observations of each species was then converted to percentages of total observations for each subestuary to normalize for differences in the frequency with which the different subestuaries were surveyed.
The SAV communities were characterized using nonmetric multidimensional scaling (NMDS). The analysis was performed in R using the meta-MDS function in the Vegan package, which performs an iterative analysis at random starts to prevent selecting a local optimum fit rather than the global optimum (Oksanen et al. 2014). We performed n = 300 iterations and limited the solution to two axes for ease of interpretation and subsequent analyses. The subestuaries were then classified into community categories using an unweighted pair-group method using arithmetic averages (UPGMA) average linkage clustering method in the Vegan Library.

Factors Affecting Community Composition
To evaluate all available covariates that might influence community composition, we fit a random forest model to each ordination axis. Random forest models are uniquely suited to dealing with large sets of predictors that may be correlated, and they can handle non-linear relationships between predictors and response variables (Cutler et al. 2007). The predictors considered were land cover (% developed, % cropland, % forest, and % wetland), shoreline condition (% bulkhead, % riprap, % developed, boat ramp density, and dock density), and subestuary characteristics (mean water depth, fractal dimension of shoreline, volume, shoreline length, estuary mouth width, estuary area, watershed/estuary area ratio, mean salinity, rate of sea level rise, tidal range, mean annual precipitation) (Patrick et al. 2014). Each random forest consisted of 500 trees constructed using a randomly selected third of the predictors and two/thirds of the observations in each tree (Cutler et al. 2007). Top predictors identified by the model were evaluated using partial dependence plots which graphically summarize the average relationship between each predictor and response variable across all other values of the covariates (Biau 2012).

Temporal Patterns of SAV Abundance and Correlates
SAV abundance through time was evaluated for the presence and slope of a monotonic trend (ΔSAV) in each subestuary using non-parametric methods because of the temporal autocorrelation in the data. We used a Mann Kendall test for the presence of monotonic trend, and slope was estimated using Theil-Sen's slope estimator from the rkt library in R (Mann 1945;Marchetto et al. 2013). To assess interannual variability in SAV abundance as a secondary response variable, we removed monotonic trends from the time series and calculated the standard deviation (SAV sd ) and the absolute value of the range (SAV ar ) for each detrended time series. Slopes were also measured for the temporal change in proportion of cropland and total development in the watershed of each subestuary from 1984 to 2006. We then evaluated our hypotheses that the % cropland, % total development, % riprap, and % bulkhead are negatively related to ΔSAV by implementing a backward stepwise multiple linear regression with 2001 land use values using the lm() function in R and removing nonsignificant predictors. We also related ΔSAV to the temporal change in land cover using a series of simple linear regressions. We used ANOVA to partition the variance in ΔSAV among the categorical explanatory variables: land cover class and community identity. To evaluate the relative importance of all potentially important covariates (land use, shoreline condition, subestuary and watershed geometry, salinity, and community composition) for which we had data, we fit another random forest model using all the covariates to predict ΔSAV. This analysis used the same set of predictor variables as the previous random forest plus the % of all observations in a subestuary that were of a given species. This same series analyses were also performed on SAV sd and SAV ar .

Temporal Changes in Community Composition and Correlates
We developed a second analysis of the ground survey observations to test for temporal changes in the SAV community in each subestuary. We limited the analysis to subestuaries that were well sampled throughout the study period by excluding subestuaries that did not have at least 10 years of data with some observations in each third (beginning, middle, and end) of the time series. Three additional subestuaries were removed from the analysis based on best professional judgment because of outlier observations and heavily uneven temporal sampling. This left 53 subestuaries for the analysis. Temporal variability was smoothed using a 10-year moving window that summarized general changes in SAV community composition through time in each subestuary (e.g., 1984-1993, 1985-1994, etc.). The smoothing approach reduced the influence of individual years and eliminated gaps in the time series. A 10-year window was used in a previous analysis of 1984-1994 data (Moore et al. 2000). For each time step, observations for the subestuaries over the 10-year period were summed and then recorded as % observations for each species for each system. This generated 19 separate subestuary by species matrices, 1 for each of the 19 full, 10-year time windows occurring between 1984 and 2012.
We adapted the previous ordination (see BCommunity Classification^section above) to assign a community to each subestuary in each year. Using the first two axes of the NMDS ordination, we mapped the ordination space occupied by each community as the area within a polygon delineated by the exterior points (the convex hull) representing subestuaries with that community. Next, we developed models to place any set of SAV species observations into the ordination space. For each of the first two NMDS axes, we fit a multiple linear regression model that estimates position along the NMDS axis using the SAV species observations as predictors. The models for the two axes were both highly significant (P < 0.0001) and explained almost all the variation in NMDS scores (axis 1: R 2 = 0.995, axis 2: R 2 = 0.982). Finally, we applied the regression models to place the species observations from each subestuary in each year into the ordination space. When a point fell within one of the community polygons in the ordination space, then the community was assigned to that subestuary in that year. Points that fell outside of the defined community cluster regions were not assigned to a cluster.
We divided the 53 subestuaries into two groups: those that maintained the same community throughout the study and those that changed from one defined community to another defined community one or more times. To evaluate if stressors affect the probability of community change, we fit a series of logistic models using land use (development and agriculture) and shoreline stressors (bulkhead and riprap) as predictors of the probability that a subestuary maintains the same community cluster through all time steps. Models were fit using a binomial distribution with the glm() function in R.

Invasive Species
To quantify the effect of stressors on invasive species, we related land use stressors (development, agriculture, riprap, and bulkhead) to the proportional composition within subestuaries of three common exotic SAV species (H. verticillata, Najas minor, and M. spicatum). The salinity range of each species was set as the observed range of salinities in which it occurred, and only subestuaries within that observed salinity range were included in the analysis for that species. Proportions of exotic taxa never approached the upper bound of 1, so each relationship was fit as a simple linear regression using the lm() function in R.

Community Classification
The NMDS ordination of the subestuary by species matrix arrived at a stable two-axis solution after 26 iterations with a stress value of 0.09 indicating a good fit for the data (Figs. 2 and 3). The hierarchical cluster analysis of the community data identified seven different communities, two of which were minor and found in only one or two subestuaries (Figs. 2 and 3). We named the five major communities we refer to as mixed freshwater (MXFR), Eurasian water milfoil (M. spicatum) and horned pondweed (Zanichellia palustris) (MZP), horned pondweed (ZP), mixed mesohaline (MXM), and widgeon grass (R. maritima) and eelgrass (Z. marina) (RZ) (Fig. 3).

Factors Affecting Community Composition
The random forest models for relating community ordination axis scores to environmental variables explained 83 and 55% of the variation, respectively, in the first and second axis of ordination. The top predictors for the first axis scores were mean salinity, rate of sea level rise, watershed/estuary area, watershed wetland proportion, and subestuary perimeter length (Fig. 4a-c). Mean salinity was clearly the best predictor, and it was positively related with the first ordination axis (Fig. 4a). The top predictors for the second ordination scores were mean tidal range, % shoreline with riprap, % developed shoreline, watershed forest coverage, and mean annual precipitation ( Fig. 4d-f). Among these, factors associated with shoreline alteration were positively related with the second ordination axis.

Temporal Patterns of SAV Abundance and Correlates
The trend analysis from 1984 to 2009 revealed differences among subestuaries. SAV abundance changed at a mean rate across subestuaries of 0.20 ± 0.46 SD (abundance year −1 ), but of the 95 systems, 36 subestuaries had significantly increasing abundance, 7 subestuaries had significant decreases, and the remaining 52 had no significant trend (Fig. 5). Across subestuaries, the mean SAV sd (standard deviation in SAV) was 3.23 and SAV ar (absolute range in SAV) was 13.18. The distributions for SAV sd and SAV ar were strongly skewed right.
Watershed development and % shoreline with bulkhead were not significantly related to ΔSAV, so they were dropped from the multiple regression model predicting ΔSAV from environmental variables. The final model predicted ΔSAV from % watershed cropland and % shoreline with riprap (Table 1). Both predictors were negatively related to ΔSAV, and the model explained 13.5% of the variation in ΔSAV among subestuaries (Table 1). The rates of change in the amount of developed land and cropland with each watershed were both positively related to ΔSAV in subestuaries (Table 2). In the two-way ANOVA model, both dominant land category and community were significantly related with ΔSAV; however, there was no significant interaction between these variables (Fig. 6,  Table 3). ΔSAV was significantly higher in MXFR communities than in any other community (Fig. 6). ΔSAV was significantly higher in forested watersheds than in human impacted watersheds (Fig. 6).
The regression models predicting SAV sd and SAV ar from environmental variables were both reduced to intercept only models because none of the predictors (cropland, development, riprap, and bulkhead) were significantly related to the response. Similarly, we observed no significant relationships between change in cropland or developed land over time and SAV sd or SAV ar (Table 1). In the two-way ANOVA models, community was significantly related to both SAV sd and SAV ar , but dominant land use type was not significant nor was there an interaction between land use and community (Fig. 4, Table 3). Mixed freshwater communities had significantly higher SAV sd than all other communities except mixed mesohaline communities and significantly higher SAV ar than all other communities (Fig. 6).
The random forest model using all covariates explained 45% of the variation in ΔSAV among subestuaries. The top five predictors in the random forest model were the proportion of the SAV community that was H. verticillata, the mean salinity of the subestuary, the subestuary's score on the first ordination axis, and the proportions of the community that were N. minor or Ceratophyllum demersum (Fig. 7a). All of these measures of freshwater or the presence of freshwater SAV species were positively related to higher ΔSAV.
The random forest models using all covariates to predict interannual variation in SAV explained 28% of the variation in SAV sd and 42% of the variation in SAV ar . The top five predictors in the SAV sd model were the proportion of the SAV community that was H. verticillata, the perimeter length of the subestuary, the amount of precipitation the estuary received, the proportion of the community that was Heteranthia dubia, and the amount of shallow water habitat (Fig. 7b). The top five predictors of SAV ar were the same, except that precipitation which was replaced with N. minor. H. verticillata, H. dubia, N. minor, and the amount of shallow water habitat were all positively related to increased interannual variation in SAV.

Temporal Changes in Community Composition and Correlates
Of the 53 subestuaries that met the data requirements for inclusion in the analysis of temporal changes in community identity through time, 42 (79%) kept the same community identity for all 19 time steps, whereas 11 (21%) changed community identity one or more times (Fig. 8, Table 4). Five subestuaries moved between the RZ community and either the MXM or ZP community. Four subestuaries moved between the MXFR community and one of the other types (Table 4). The probability of changing SAV community identity was significantly higher in subestuaries with more riprapped shoreline but was not significantly related to the proportions of cropland or developed land in the local watershed of shoreline bulkhead (Table 5).

Invasive Species
The observations of exotic species in the ground survey data indicate that subestuaries with salinities below 5.5 support H. verticillata and N. minor, while M. spicatum occurs in subestuaries up to 9. Among subestuaries with 0-5.5 salinity, watershed agriculture, and % riprap shoreline had positive relationships with M. spicatum and negative relationships with N. minor, while H. verticillata was not significantly related to any of the predictors (Table 6). M. spicatum was not related to human stressors among subestuaries with 5.5-9 salinity (Table 6).

Discussion
We characterized the SAV communities in Chesapeake Bay subestuaries and analyzed the relationships among community, terrestrial stressors, and changes in SAV abundance and community identity over time. To our knowledge, this is the first effort to integrate all these aspects of estuarine SAV dynamics into a single set of analyses. Salinity and terrestrial stressors strongly influenced SAV community composition (Fig. 4), which, in turn, affected trends and interannual variability in SAV abundance (Fig. 5) and the impacts of stressors on temporal variability (Tables 1, 2, and 3). As we hypothesized, human land use (the proportion of cropland in the local watershed and the amount of riprap shoreline) had a negative effect on the temporal SAV trend in subestuaries (Table 1). Surprisingly, developed land did not have the hypothesized negative effect on SAV trends (Table 2). While subestuaries composed of species tolerant to stressors like low light availability had the overall largest positive trends; these same subestuaries also exhibited the strongest negative responses to terrestrial stressors. Additionally, we did not find the expected relationships between the stressor intensity and changes in community or the prevalence of exotic SAV species (Table 6). Interestingly, the presence of exotic H. verticillata was the single best predictor of SAV temporal trend.

SAV Communities
Our community classifications are similar to those reported by Moore et al. (2000). Our mixed mesohaline community is comparable to their Potomageton community, and both analyses identified a mixed freshwater and a Z. marina/R. maritima community. However, there are some key differences in the classifications. Moore et al. (2000) identified a R. maritimedominated community that we did not observe, and we identified an M. spicatum/Z. palustris as well as a Z. palustris community that Moore did not observe. The primary reason for the latter difference is that Moore et al. (2000) considered the entire bay including mainstem river and bay habitats where M. spicatum and Z. palustris typically co-occur, while we analyzed communities only within subestuaries (Fig. 1), including some subestuaries where Z. palustris occurs by itself. The reason that subestuaries may differ from mainstem is that while Z. palustris tolerates salinities as high as 10-15 (Batiuk 2000), Z. palustris is unlikely to appear in the Chesapeake Bay mainstem or mainstem rivers because of its shallow root system. Our analysis also included 17 years of additional ground observation data (1996)(1997)(1998)(1999)(2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012), and differences between the studies may reflect actual changes in species distributions. We did observe that nearly a quarter of the subestuaries in the analysis changed from one community cluster to another at least once between 1984 and 2012 (Table 4).
Salinity tolerances largely explain how species are arranged along the first ordination axis, but the second axis appears to represent anthropogenic factors because riprap, watershed forest, and bank development were three variables most associated with that axis (Fig. 4). Most of the variation along the second axis was in the middle and lower salinity zones (Fig. 4), which is consistent with the higher species diversity there than in the high salinity zone (Moore et al. 2000). The mesohaline zone has four different community clusters, and there are large differences among them in local watershed land cover and shoreline land use. For example, there is a higher proportion of watershed development and shoreline bulkhead in the M. spicatum/Z. palustris cluster (Fig. 3) than in the other mesohaline clusters. This may be partially driven by M. spicatum's positive association with % riprap within the 0-5.5 salinity range (Table 6). Interestingly, Z. palustris exists in monoculture for a cluster of subestuaries. Fig. 4 Results of random forest models relating the first two ordination axes to potential explanatory variables. a, d Variable importance plot for the top 5 predictors of NMS1 and NMS2, respectively. b-f Relationships between the two ordination axes and the top 2 predictors of each axis, shown as partial dependence plots (left) and bivariate plots of the raw data (right) Changes in community identity occurred only in the upper and middle bay, not in the polyhaline zone. Z. marina and R. maritima are the only two Chesapeake Bay species that can survive and grow at high salinity levels, so polyhaline subestuaries cannot change to other communities containing different native species (Batiuk 2000). However, more than half of the observed community changes outside the polyhaline were transitions to and from the R. maritima/ Z. marina community in the mesohaline zone. R. maritima exhibits poorly understood boom and bust cycles (VIMS ground survey), which contribute to the community changes in the mesohaline zone. Fig. 5 The spatial patterns of the temporal trend ΔSAV (a) and temporal variation SAV ar (b) of each subestuary. SAV ar is the absolute range of the detrended variation in SAV abundance observed in each system. SAV sd is not shown because it has a nearly identical pattern to SAV ar

Temporal Variability in SAV Abundance
The average trend among subestuaries is positive, indicating that SAV is recovering, as reported in previous bay-wide analyses (Orth et al. 2010). However, the variation in trend among subestuaries is high, and two thirds of the subestuaries had no trend or a declining trend. Among subestuaries, cropland and riprapped shoreline were negatively associated with SAV recovery but had less effect than community. This supports the hypothesis that stressor impacts on SAV differ among communities (Patrick and Weller 2015). Unexpectedly, the negative impacts of human stressors on trends were strongest in the mixed freshwater community, counter to our hypothesis that freshwater SAV communities would be less sensitive to human impacts. The mixed freshwater community also had higher interannual Fig. 6 Relationships among temporal measures of subestuary SAV abundance (ΔSAV, SAV sd , and SAV ar ), community, and dominant land use in the local watershed. Human land use includes systems dominated by either cropland or developed land fluctuations than any other community except the mixed mesohaline community (Fig. 6), matching similar observations in Gulf of Mexico estuaries (Martin and Valentine 2012). Freshwater and mesohaline communities, particularly those with R. maritima and H. verticillata, can exhibit boom and bust cycles (Cho et al. 2009;McChesny 2010). Across all the Chesapeake Bay subestuaries, SAV recovery rates are highest for the mixed freshwater communities that occur in subestuaries with natural or mixed land use watersheds (Fig. 6). Previous research has documented that after SAV beds reach a critical threshold in patch size, recovery can accelerate through positive feedbacks ( Similar positive feedbacks have not been observed in the mesohaline or polyhaline parts of Chesapeake Bay. The documented potential for a critical size effect in mixed freshwater communities may explain why SAV recovery is stronger for mixed freshwater communities than in the other communities. Positive feedbacks could be overcoming the negative effects of human land use in freshwater subestuaries. The high rates of recovery in mixed freshwater communities also explain why the rate of change in the proportion of developed land is positively related to the rate of change of SAV abundance. The subestuaries with the highest rates of watershed development were freshwater systems with low mean development coverage. Subestuaries with watersheds that are already heavily developed already have severely reduced SAV and are not changing as rapidly. Fig. 7 Results of a random forest models relating temporal patterns to potential predictors. Response variables are the change in SAVabundance through time (a ΔSAV) and the standard deviation of detrended SAV abundance through time (b SAV sd ). Each column (a-b) has the variable importance plot for the top 5 predictors at the top (1), followed by graphical depictions of the relationships between temporal patterns and the top predictors (2-6), shown as partial dependence plots (left) and bivariate plots of the raw data (right)

Factors Affecting Community Fidelity
One fifth of the subestuaries changed community identity at least once. There was a significant negative relationship between the probability of community change and the amount of riprapped shoreline. Riprap is predicted to decrease habitat suitability by increasing depth and lateral scour (Prosser et al., this volume). Loss of shallow water habitat and strong scour effects could lead to frequent local extirpation and colonization cycles with community changing due to change invasions and priority effects (Alford and Wilbur 1985;Chase 2003). In communities where multiple stable communities are possible (i.e., multiple stable equilibria), the identity of the first colonists exerts a strong influence on the identity of the community that develops. We also observed that there were no stable communities that occurred in the 3-11 salinity range. Many species are near the upper or lower bounds of their salinity tolerance in this range, and slight interannual variation in salinity could lead to either loss or increases in abundance (Batiuk 2000;Kemp et al. 2004).
Unique Behavior of Z. palustris Z. palustris differs from other SAV in the Chesapeake Bay in two ways: It is strictly an annual species, and it completes its life history between fall and spring. Seeds germinate in the fall; plants grow rapidly in the spring and then flower and die by early July before the aerial surveys so that subestuaries dominated by Z. palustris group appear devoid of SAV in the aerial surveys (Davis 1985;Lombardi et al. 1996). Paleoecological research has documented that subestuaries of Chesapeake Bay (such as the Middle River, Back River, and Rock Creek off the Patuxent River) had large increases in Z. palustris in the 1700s and 1800s following European settlement and land clearing (Brush and Hilgartner 2000). However, the question remains, why is Z. palustris successful in these systems while other species are not?
Dispersal limitation could help explain the presence of Z. palustris in systems lacking other freshwater species Some subestuaries always remain the same community, while other subestuaries switch communities one or more times (Kendrick et al. 2012;Lloyd et al. 2016). Z. palustris does well in freshwater streams worldwide (Wilby et al. 1998;Stegen et al. 2000), so reservoir populations in the streams flowing into a subestuary could support subestuary Z. palustris. For example, the Rhode River subestuary is devoid of SAV in the aerial surveys but has an annually recurring population of Z. palustris in the upstream Muddy Creek which is the major tributary of the larger subestuary basin (personal observation C. Patrick). Those reservoir populations could provide a steady supply of propagules, whereas other species without reservoir populations have been dispersal limited and have not re-established despite adequate water quality conditions. An alternative explanation for the absence of other species in the subestuaries supporting only Z. palustris is that water clarity in winter and spring is high enough to support plant growth, but not in summer when other SAV species need to complete their live cycles. Furthermore, Z. palustris can thrive under eutrophic conditions that other species may find stressful (Van Viersson 1982). A final contributing explanation for pattern is that Z. palustris has a wide salinity tolerance (Batiuk 2000), so it may persist where salinity oscillates so widely that both freshwater species and higher salinity species are excluded.

Invasive Species
Exotic SAV species have been in the Chesapeake Bay for decades, and their effects on community dynamics and SAVabundance are an area of active research (Orth and Moore 1984;Posey et al. 1993;Rybicki and Landwehr 2007). Where restoration goals are far from being met, it may be tempting to suggest that exotic SAVs are more desirable than no SAV at all. The SAV restoration goals of the CBP focus on the area restored (Batiuk 2000;Orth and Wilcox 2009), not on its species composition, and the aerial survey used to monitor SAV does not distinguish species, so gains in SAV coverage from exotic taxa are implicitly counted as positive. Two of the common exotic species in Chesapeake Bay, M. spicatum and H. verticillata, can increase the SAV abundance in freshwater systems to nuisance levels that obstruct boat traffic (Langeland 1996). There was concern that M. spicatum would overrun Chesapeake Bay's freshwater and low salinity SAV habitat, and it did undergo a period of the massive expansion in the 1950s. It subsequently died back in most areas and has not regained its prior abundance levels (Orth and Moore 1984).
Our hypothesis that the presence of exotic species would be related to higher variability in SAV abundance was not supported, and different exotic species had different effects of SAV abundance. M. spicatum abundance was positively related to the cropland and riprap shoreline among the 0-5.5 subestuaries ( Table 5), suggesting that it is stimulated by anthropogenic disturbances, as reported in lakes (Trebitz and Taylor 2007). However, while M. spicatum tends to dominate other macrophytes in lakes, it is not outcompeting the other taxa in the fresher waters of Chesapeake Bay. M. spicatum is outcompeted by Vallisneria americana in high-energy zones in the Gulf of Mexico (Martin and Valentine 2012), suggesting that M. spicatum is not well adapted to high-energy open water systems. N. minor, another exotic species that can thrive in impacted environments (Trebitz and Taylor 2007), was also negatively related to cropland and riprap shoreline (Table 5), but its presence was significantly related to increases in SAV abundance over time. The dependent variable is the proportion of all SAV observations within a subestuary that were the invasive species. Only subestuaries within the indicated salinity ranges were included H. verticillata was strongly related to SAV recovery rate (Fig. 7), and the proportion of H. verticillata explained 53% of the variation in SAV temporal trend among all subestuaries. However, the abundance of H. verticillata was not related to terrestrial stressors (Table 6). H. verticillata never made up more than 37% of the total ground observations for any one subestuary and average 19 ± 13% SD of the observations in subestuaries where it was present. These percentages are similar to those reported from more quantitative community surveys in the Potomac River (Rybicki and Landwehr 2007). Finally, more dominant species like V. americana and the Potamogeton species decline as H. verticillata increases, but rarer species like C. demersum, H. dubia, Najas guadalupensis, and Najas gracillima become more common. H. verticillata may be acting as a pioneer or foundation species, stabilizing sediments and facilitating the establishment and spread of some of the less common SAV species in Chesapeake Bay (Rybicki and Landwehr 2007).
Experimental work has demonstrated that native V. americana can restrict the colonization of H. verticillata if the native community can reduce water column nutrient concentrations enough to prevent H. verticillata establishment (Chadwell and Engelhardt 2008). The high nutrient requirements of H. verticillata may preclude a comparable exclusion effect on native taxa when H. verticillata establishes first. Unpublished experimental work also suggests that H. verticillata cannot outcompete native V. americana (McChesny 2010). Manipulative restoration experiments planting mixed beds with and without H. verticillata would help clarify whether the observed associations with H. verticillata are correlative or causal.

Management Implications
Our study enhances overall understanding of the factors contributing to spatial variation in SAV recovery and that knowledge can improve management strategies for SAV restoration and conservation. Our results indicate than human activities have inhibited the rate of SAV recovery or even prevented SAV recovery in some subestuaries. However, the effects of stressors on SAV recovery and variability differ strongly among communities, supporting the idea of customizing management strategies to community composition (Kemp et al. 2004;Patrick and Weller 2015). Some general lessons for SAV management also emerged. In middle to low salinity subestuaries, the amount of riprap shoreline was a predictor of differences in SAV community among subestuaries and of the probability of community change in a subestuary. Riprap shoreline was also a negative predictor of SAV recovery rate. These findings add to the growing body of literature on the negative effects of shoreline alteration on nearshore macrophytes and seagrasses (Strayer et al. 2012;Patrick et al. 2014;Patrick et al. 2016). We interpret the results to mean that management efforts to slow the rate of shoreline armoring or to replace hard armoring with marsh and Bliving shorelines^are likely to have positive impacts on SAV.
We also suggest that the dynamics of particular SAV species merit further investigation. For example, H. verticillata is widely considered an undesirable species but was strongly related to SAV recovery in our analysis and does provide important ecosystem services (Moxley and Langford 1982;Posey et al. 1993;Rybicki and Landwehr 2007). More research is needed to understand the role that H. verticillata plays in SAV communities. Hydrilla was first recorded in the Potomac River in the early 1980s, has spread rapidly, and now is persistent in SAV communities throughout low salinity regions of the bay. Similarly, the occurrence of Z. palustris in subestuaries that lack other species and have low SAV abundance is interesting and merits further investigation in the field. Better understanding of the interspecific interactions and the dynamic changes in community composition and structure over time may lead to improved understanding of how to achieve higher restoration planting success.
Our findings support the general idea that salinity is the first control on SAV community identity; however, within the physiological bounds of salinity zone, community composition can change in response to other environmental drivers, which can alter the temporal pattern of SAV abundance in a subestuary. A greater appreciation for the naturally dynamic nature of community composition and its role in regulating SAV abundance dynamics will improve our conservation and management of these fundamental estuarine ecosystems.