The IRL is a shallow (mean depth < 2 m), bar-built lagoon that stretches over 250 km along the east coast of Florida from Ponce de Leon Inlet in the north to Jupiter Inlet in the south (Figure 1). The system’s biological diversity can be explained primarily by the fact that it spans the boundary between the temperate Carolinian Province and the subtropical Caribbean Province (Gilmore, 1995). In addition, diversity is affected by variation in hydrodynamic characteristics caused by the length of the system and the locations of five inlets and one navigational lock. Ponce de Leon Inlet is at the extreme northern end of Mosquito Lagoon; the lock at Port Canaveral intermittently connects the ocean to the Banana River Lagoon; the next connection is 140 km to the south at Sebastian Inlet; and Fort Pierce, St. Lucie, and Jupiter inlets are all found in the remaining portion of the lagoon. The three southern inlets increase oceanic exchange in that region (Smith, 1993). In contrast, the lack of inlets in the northern portion of the system leads to residence times that can exceed 1 year and weak tidal forcing, with tidal amplitudes being less than 10 cm. In fact, seasonal winds and oceanic forcing driven by changes in the Gulf Stream cause larger variations by raising water levels 5–30 cm throughout the IRL (Smith, 1993, 2001; Woodward-Clyde Consultants, 1994).

For this paper, the IRL is subdivided into smaller areas labeled reaches (Figure 1). The nine reaches are delineated by evaluating similarities in time series of water quality parameters (Lasi et al., this volume). Reaches 1–4 are microtidal (<10 cm), with expansive, shallow shelves suitable for seagrass (Smith, 1987). In addition, these northern reaches are adjacent to some of the most pristine areas within the IRL system, including Merritt Island National Wildlife Refuge, Canaveral National Seashore, and multiple areas designated as Outstanding Florida Waters. Throughout the system, the bathymetry and relatively stable and favorable water quality have supported persistent coverage of seagrass.

Seagrass is sampled at two scales: large-scale, lagoon-wide mapping from aerial imagery and finer-scale, fixed transects. Thirteen maps spanning the whole lagoon document conditions in 1943, 1992, 1994, 1996, 1999, and every 2 years from 2003 to 2019. Surveys of fixed transects are used to detect small-scale changes within the lagoon’s seagrass beds, and they began in 1994.

Mapping is based on visual interpretation of aerial images, primarily at a 1:24,000 scale but in some cases at a 1:10,000 scale. Features on the aerial images are identified with the aid of stereoscopic analysis, photo-interpretation keys, and ground truthing. Features are classified according to slightly modified Florida Land Use, Cover and Forms Classification System codes (Florida Department of Transportation, 1999). The resulting polygons are connected to create a layer showing the areal extent of seagrass. The accuracy of classifications has been evaluated from 1999 onward, and it averages 80% or higher.

Starting in 1994, 73–97 fixed transects placed throughout the lagoon’s seagrass beds have been surveyed at least twice a year, during summer and winter, to approximate times of maximum and minimum abundance within years (Virnstein and Morris, 1996). Each transect is delineated by a graduated line extending perpendicular to the shore out to the deep end of the seagrass canopy (EOC). For consistency with mapping, the EOC is designated according to the distribution of larger species that form the canopy (i.e., Halodule wrightii, Syringodium filiforme, Thalassia testudinum, and Ruppia maritima) because the three Halophila species found in the lagoon are small, ephemeral, and difficult to map accurately from aerial images. The EOC is defined as the point beyond which percent cover in 1-m² quadrats drops to zero or remains < 1% for 30 m. At predesignated points along transects, a suite of standardized, non-destructive measurements are made within 1-m² quadrats divided by strings into 100 cells, each 10 cm by 10 cm. Pertinent measurements are: (1) visual estimates of percent cover of all seagrass, (2) transect lengths or distance from shore to the EOC, (3) depths where quadrats are sampled, and (4) depth at the EOC (Morris et al., 2001).

Light is mandatory for growth and survival of seagrasses, and light availability is affected strongly by depth. Depths in the IRL are influenced by seasonal changes and by tides when in proximity to the inlets. For this study, depths measured in the field are adjusted to depths below mean water level (MWL) using depth soundings taken throughout the IRL at 15.2-m intervals along east-west tracks that were spaced approximately 150–300 m apart (Coastal Planning and Engineering, 1997). The resulting bathymetry is based on over 230,000 depths that are referenced to North American Vertical Datum of 1988 and adjusted to depth below MWL.

The depths below MWL yield isobaths that are used to adjust depths measured in the field. Based on a linear interpolation between measured depths, the depth at the intersection of the transect and the 0.75-m isobath is calculated for each transect in each survey. The appropriate difference between the predicted depth and 0.75 m is used to correct all relevant field measurements. The resulting depths below MWL are rounded to the nearest 0.1 m.

The distributions of seagrasses are influenced by water temperature, salinity, and attenuation of photosynthetically active radiation (PAR, wavelengths of 400–700 nm). To assess these influences, relevant data are extracted from monthly sampling at fixed stations near the transects.

Water temperatures and salinities are measured in the field with multiparameter instruments equipped with the relevant probes. The instruments are calibrated and verified using standards and procedures that agree with state guidelines.

Since 1996, the coefficient of attenuation for PAR (K_d) is determined by applying the Lambert-Beer equation to simultaneous readings recorded by three spherical quantum sensors. Two of the sensors are held at 0.2 and 0.5 m below the water’s surface, and the third is held 0.3 m from the bottom in water ≤ 1.8 m deep or at 1.5 m from the surface in deeper water (Morris et al., 2002). For the purpose of this paper, the Lambert-Beer equation is rearranged to yield the percent of subsurface irradiance reaching the depth of interest.

Halodule wrightii is the most common and abundant seagrass species in the IRL (Morris et al., 2021), so the stores of carbon (C), nitrogen (N), and phosphorus (P) in its aboveground and belowground tissues are evaluated. The evaluation is based on triplicate samples collected quarterly in 2014 and 2015 at mid-bed (0.5–1.0 m deep) near three transects in the northern Indian River Lagoon (reaches 2–3) and three in the Banana River Lagoon (reach 4). Processing samples involves removing epiphytes and debris, rinsing briefly to remove excess salt, partitioning aboveground and belowground tissues, and drying at 80°C (Hanisak, 2016).

Dried samples ground to fine powder with a mortar and pestle are stored in capped glass vials at room temperature. Subsamples in small vials are held at 80°C overnight to ensure they are dry before being shipped to the University of Maryland’s Nutrient Analytical Services Laboratory for determination of C, N, and P contents. Analytical procedures, quality control/quality assurance measures, and other details can be found at: https://www.umces.edu/nutrient-analytical-services-laboratory.

Areal extents, distance to the EOC, depth at the EOC, percent cover, data on water quality, and elemental contents of H. wrightii are analyzed to identify spatiotemporal changes in the quantity and quality of seagrass beds, the location and characteristics of persistent patches of seagrass, drivers of changes in seagrass, and consequences for the cycling of carbon, nitrogen, and phosphorus. In some analyses, early years are excluded and reaches 7, 8, and 9 are combined or excluded because certain data are lacking.

Areal extents in hectares in each of seven reaches (reaches 7, 8, and 9 are combined) for the 11 years with maps are used in a permutation analysis of variance (PERMANOVA; Anderson et al., 2008). The analysis is based on Bray-Curtis resemblance measures calculated with range standardized data. The model treats reach as a fixed factor and years of mapping as repeated measures nested within reaches. Furthermore, a hierarchical cluster analysis based on Bray-Curtis resemblance measures and single linkage clustering is used to identify years with similar areal extents across all reaches.

In an effort to evaluate the extent of seagrass in years without maps, Pearson correlation coefficients are calculated for combinations of range standardized areal extents, distances to EOCs, and standardized depths at the EOCs. Strong, positive correlations provide increased confidence in using measures from transects to evaluate changes in seagrass extent between maps.

Persistence is defined as the number of consecutive maps in which a given polygon of seagrass is present. Persistence is determined for each reach across the 11 years with maps between 1996 and 2019. The polygons are intersected with bathymetry to identify the depths relative to MWL spanned by each polygon in each combination of reach and score for persistence (Figure 2). Depth zones with more persistent seagrass are determined from summary statistics. In addition, mean percent covers from 1996 to 2019 are calculated for 0.1-m depth increments for transects in reaches 1–6, with the resulting data analyzed with a hierarchical cluster analysis based on Pearson correlation similarity measures and group average clustering (Anderson et al., 2008).

Percent cover values from summer surveys of transects in reaches 1–6 are evaluated in two ways. The mean percent cover recorded in quadrats placed within the beds is calculated for the whole system as a measure of the quality of seagrass within the beds. The same data are adjusted to account for changes in distances to the EOCs so that they accurately reflect the amount of seagrass in the system, i.e., a small area with dense seagrass is equivalent to a larger area with sparse seagrass. The adjustment for each transect involves adding values of zero beyond the recorded EOC out to the maximum distance recorded from 1996 to 2019. The resulting data are range standardized and used in a PERMANOVA based on Bray-Curtis resemblance measures. The model treats years and reaches as fixed factors, depth zones as a random factor nested in combinations of years and reaches, and transects as repeated measures nested within combinations of years and reaches.

Changes in the distribution and percent cover of seagrass are driven by changes in salinity, water temperature, and light availability. Environmental limits applied here represent mean values for H. wrightii, S. filiforme, and T. testudinum, i.e., stress as indicated by reduced mean percent cover is caused by salinities below 23 psu, water temperatures below 20°C, and percent of subsurface irradiance below 27% at the deep edge of the range of depths where seagrass persisted longest (Morris et al., 2021). These limits reflect changes in production or loss of leaves in response to salinity, temperature, and light availability, but they also reflect responses to interactions among these factors and the results of intraspecific and interspecific interactions. Using water quality data collected from 1997 to 2018, mean values for each parameter are calculated for each combination of reach and month within year. Heatmaps identify reaches and months where seagrasses experienced stress. In addition, mean percent cover values from surveys are calculated for each transect. The difference in cover between each pair of consecutive surveys is calculated. Those differences are culled to a set that exhibited decreases in cover following at least 1 month of subsurface irradiance below 27%, and the resulting data are converted to proportional changes in cover per month of low light before being ordered by increasing initial mean percent cover. A boundary analysis using squared Euclidean distances calculated for moving windows encompassing four, six, and eight differences identifies a threshold in standardized proportional differences and an associated initial mean percent cover (Ludwig and Tongway, 1995). The validity of the threshold is tested with a t-test that compares the proportional differences above the threshold to those below it without assuming equal variances. A significant t-test indicates that the initial mean percent cover associated with the threshold represents an ecological boundary in resilience to low light.

Changes in the areal extent or percent cover of seagrass affect the cycling of C, N, and P. Given relatively constant loads of these elements to the IRL, decreases in the amount of seagrass due to mortality or decreases in uptake due to reduced growth result in more C, N, and P being available to other primary producers, including phytoplankton and drift macroalgae. To elucidate such changes in availability for these key elements, the relevant percent cover values, as determined by statistically significant terms in the PERMANOVA described above, are converted to biomass (g dry weight m^–2) using Equation 1 (Hanisak, 2001), and the resulting biomasses are converted to amounts of C, N, and P using the percentages of each element found in the aboveground portions of H. wrightii shoots. Prior to estimating the effects on cycling of C, N, and P, temporal variation in the content of aboveground and belowground tissues is investigated via a series of PERMANOVAs that treat years as a fixed factor and sampling events as a random factor nested in years.

Across all reaches, approximately 3,000 ha of seagrass were lost between 1943 and 1994 (Figure 3 and Supplementary Table 1). The timing of this loss cannot be determined rigorously due to a lack of data. However, Haddad and Harris (1985) noted approximately a 35% decline in seagrass in the 1950s that correlated with decreased water transparency, and dredge and fill projects removed seagrass to create new residential and recreational areas throughout the lagoon prior to the 1970s (Woodward-Clyde Consultants, 1994). Therefore, the expansion of seagrass from the 1990s to 2007–2009 represented a recovery. As the extent of seagrass beds expanded from 2003 to 2009, the mean percent cover within the beds decreased, whereas the mean percent cover standardized to maximum transect lengths remained relatively stable (Figure 3). These results indicated that a similar biomass of seagrass was distributed more sparsely as the beds expanded.

A PERMANOVA indicated significant differences in the areal extent of seagrass among the years with maps after 1994 (Table 1). A hierarchical cluster analysis indicated three groups of years with similarities of 91% (Figure 4). The years from 1994 to 2009 formed one group; 2011, 2013, and 2015 formed a second group; and 2017 and 2019 formed a third group. These groups tracked the lagoon-wide increase in hectares of seagrass from 1994 to 2009, the initial decrease in areal extent in 2011, and the second decrease in extent in 2017 (Figure 3).

In most reaches, correlations indicated useful relationships among distance to the EOC or transect length, depth below MWL at the EOC, and areal extent (r ≥ 0.70; Table 2). The lack of data on depths precluded two correlations in reaches 7–9, but those reaches, along with reach 1, had the weakest correlations.

Overall, distance to the EOC and depth at the EOC paralleled changes in areal extent. As an example, the consistent mean distances to the EOC suggested there was little change in areal extent in years between maps and relatively large decreases in areal extent in reaches 1–6 and 8 during 2010–2011 and 2016–2017 (Figure 5).

In combination with bathymetry, maps also provided a means of assessing persistence of seagrass at different depths. Box plots of the depths intersected by polygons that persisted across differing numbers of consecutive maps from 1996 to 2019 indicated that polygons detected in seven or more maps tended to span depths of 0.5–1.0 m (Figure 6).

Persistence varied among reaches (Figure 7). Reaches 1 and 2 had the highest percentage of hectares that persisted across 10–11 maps (66.7 and 37.8%, respectively). In reach 3, no hectares persisted for more than seven consecutive mapping years, which meant seagrasses did not persist beyond 2011. Reaches 4 and 5 had the fewest persistent hectares, with only 5% of them present for 10 or more years, and reaches 6 and 7–9 had a moderate number of persistent hectares, i.e., 15–25%.

In addition to varying in extent and location, seagrass beds also varied in quality, which was documented as percent cover of seagrass within mapped footprints. Cover varied both temporally and spatially (Figure 8). Cover was more stable in the northern reach (reach 1) and in three southern reaches (reaches 7–9 combined) until 2016. In contrast, the central reaches, 2–6 were more variable and exhibited larger losses in mean percent cover after 2011. Interestingly, percent cover decreased following three hurricanes in 2004–2005, but it rebounded in most reaches.

A boundary analysis of proportional decreases in mean percent cover for transects standardized to the number of months with light below 27% of subsurface irradiance indicated that losses were greater when initial mean percent cover was below 20% (one-tailed t₁₀₂ = −16.64; p < 0.001). On average, the rate of decrease in mean percent cover per month of exposure to low light was approximately 2 × faster when the initial mean percent cover was below 20%.

Given the variation in persistence across depths, a hierarchical cluster analysis was performed on mean percent covers in different depth zones below MWL. The analysis yielded three groups: shallow (0.1–0.4 m), mid (0.5–0.9 m), and deep (1.0–2.0+ m) at 43% similarity (Figure 9).

Further evidence for significant variation in cover among the three depth zones came from a PERMANOVA that indicated mean percent cover during the summer varied significantly among depths in the combinations of year (1996–2019) and reach (reaches 1–6; Table 3). Prior to 2011, the mean percent cover at the mid depths (0.5–0.9 m) tended to be higher (Figure 10). During the period of record, all reaches exhibited declines in mean percent cover, with mean cover in reaches 3, 4, 5, and 6 declining in 2011 and mean cover in reaches 1 and 2 declining in 2016–2018 (Figure 10).

Based on the mean of the environmental limits established for H. wrightii, S. filiforme, and T. testudinum (Morris et al., 2021) and the persistence of seagrass at depths of 0.5–0.9 m, limits for stress were set at salinities below 23 psu, water temperatures below 20°C, and percent of subsurface irradiance reaching 0.9 m below 27%. Heatmaps indicated that longer periods of stressful salinities affected reaches 3–5 in 1997–2000, 2002–2003, 2004–2006, 2008–2009, and 2017–2018 (Figure 11). Temperatures were stressful during the winters of most years, with longer periods of stress in the winters of 1998, 2004, 2008, 2009, 2010, and 2011 (Figure 11). Stress from reduced light availability occurred repeatedly in reach 1, and more widespread stress occurred in 2010–2011, 2016, and 2018 (Figure 11).

Stress was evaluated as a change in cover, and it did not equate to mortality. The periods of low salinity led to decreases in cover of 0.3–2.8%. The periods of low temperatures led to 1.3–6.8% decreases in cover. In contrast, the periods of low light availability led to decreases in cover of 3.8–12.9%, with the largest decrease related to the 2010–2011 event.

Climate and weather affected salinities and temperatures, and internal conditions affected the availability of light. Low salinities corresponded to two El Niño periods with relatively high rainfall, a period with multiple hurricanes, the time when Tropical Storm Fay hit the region, and the passage of Hurricane Irma. Air temperatures were at or below the 10th percentile (8°C) for more days leading into and during periods of temperature stress. Stress from reduced light availability coincided with 3–5 fold increases in phytoplankton biomass (Phlips et al., 2021).

Given the periods identified as stressful and the changes in areal extent and cover of seagrass, it appeared that decreased light penetration was the major influence on seagrass. Areal extent and mean percent cover decreased from 2009 to 2011 and from 2015 to 2017, which bracketed periods with widespread stress due to decreased light availability at 0.9 m. In contrast, areal extent and cover increased in 2013 and 2015, when all stressors were absent.

Changes in the amount of seagrass present in the IRL would have altered the cycling of C, N, and P. Analyses of tissues from H. wrightii indicated that there were statistically significant differences in aboveground %N and %P and belowground %C, %N, and %P among sampling events, with percentages increasing during the non-growing season and decreasing thereafter (Figure 12; PERMANOVAs; p < 0.001 for significant tests and p > 0.05 for the non-significant test). Mean %C in aboveground tissue did not vary significantly. For all stations, aboveground %N was always greater than belowground %N.

Changes in C, N, and P stored in the aboveground tissues of H. wrightii were estimated by applying Equation 1 (Hanisak, 2001) to the mean percent cover for every transect in every sampling event to yield biomass, and converting those values to the mass of each element using Equations 2, 3, and 4, which were based on the overall mean percentage content of all aboveground tissues sampled. Non-destructive surveys of transects did not yield data on belowground biomass, so changes in storage of C, N, and P were underestimated if belowground tissues were lost.

The estimated elemental contents of the aboveground tissues of H. wrightii were extrapolated to the aerial extent of seagrass by multiplying the relevant value by the number of hectares in each reach, with values for years without maps based on the area from the previous year (Figure 13). The amount of material stored in seagrass beds declined in 2011 in reaches 2, 4, 5, and 6, and it declined later in the other reaches. By 2019, the amount of C, N, and P stored in seagrass had decreased by over 88% in reaches 1–5 and by over 58% in reaches 6, 7–9.