Introduction

Understanding the effects of altering the timing and volume of natural freshwater flow to estuaries is one of the most pressing challenges for coastal scientists (Alber 2002; Estevez 2002; Gillanders and Kingsford 2002; Davis et al. 2005; Erwin 2009). The competing needs of agricultural irrigation, municipal use, wetland drainage, flood control, and maintenance of ecosystem services such as aquifer recharge, biodiversity conservation, wildfire suppression, nutrient management, and healthy fisheries create a challenge for resource managers tasked with harnessing freshwater as it flows from upland landscapes to coastal ecosystems. Some needs require extraction or redirection of water resources, others require fresh water to remain in the watershed. Compounding this problem, water management decisions are now more of a moving target than ever due to the accelerated development of coastal zones and changes in climate. To address this need, it is important to understand the changes in estuarine systems over long timescales and under various flow scenarios.

Some causes of increased freshwater flow to estuaries include deforestation which promotes rapid runoff and less infiltration into groundwater, draining wetlands with trenches and canals, channelization and dredging of rivers, and water control structures that combine runoff from multiple basins into a single outfall (Grange et al. 2000; Gillanders and Kingsford 2002; Bellio and Kingsford 2013). More freshwater flow often results in a cascade of altered bio-physical responses such as more carbon, nutrients, phytoplankton, zooplankton, freshwater tolerant fishes, and other biota (Grange et al. 2000; Palmer et al. 2011; Bittar et al. 2016; Loh et al. 2017), but may also cause sedimentation and smothering, eutrophication and hypoxia, pollution, and lack of natural drying cycles necessary to maintain ecosystem function (Rehage and Loftus 2007; Bellio and Kingsford 2013).

Reduced freshwater flow to estuaries can be caused by agricultural and municipal use, impounding or blocking routes of natural flow, or diversion of watershed drainage patterns from one outfall to another (Alber 2002; Gillanders and Kingsford 2002; Lorenz 2014). Less flow often results in reduced nutrients and primary productivity (Ley et al. 1994), loss of sediment that maintains estuarine morphology, shifts in the distribution of vegetation (Ross et al. 2000), encroachment of higher salinity, and a chain of other ecosystem effects. These may include an increase in saltwater tolerant species, reduced cues for larval recruitment (Grange et al. 2000; Smith et al. 2008; Criales et al. 2010), changes in prevalence of some pathogens (Volety et al. 2014), or even changes in density of predators at the top levels of the food web (Palmer et al. 2011; Lorenz 2014).

Altered flow not only refers to the amount of fresh water reaching the estuary but also the timing of its arrival (Alber 2002; Estevez 2002; Palmer et al. 2015). Natural wet/dry seasonal cycles in freshwater delivery or the timing of extreme events such as hurricanes and drought may be just as important to estuarine communities as the amount of water delivered (Rehage and Loftus 2007; Serafy et al. 1997).

Whether there is more or less flow, effects are often estuary-specific and depend on initial estuarine morphology, trophic structure, and the relative amount and seasonal timing of natural versus altered flow (Estevez 2002). Depending on how flow is altered, community structure may change toward a suite of more salt- or freshwater-adapted taxa, or the species present may remain consistent, but their relative abundance may change to favor those better adapted to the altered flow (Alber 2002; Estevez 2002).

There are several challenges that limit forming broad generalizations from the many studies on this topic. Most studies span only a few years, which may not represent either the average or range of interannual conditions that shape biotic communities. Many studies only investigate either increased or reduced freshwater flow, seldom both, thereby limiting opportunities to fully understand environmental drivers. Still other studies in riverine estuaries sample along the upstream/downstream salinity gradient; however, this approach adds further correlated variables such as distance to the ocean and changes in geomorphology and habitat that complicate the interpretation of environmental effects. Lastly, most studies lack a control system with relatively natural flow for comparison. These limitations often occur because sampling design is constrained by the conditions in each landscape being altered, a setting that is almost always beyond the control of the experimenter. Together, these circumstances have limited the direction (e.g., more or less fresh water) and magnitude of altered environmental conditions, increased potential bias caused by a period of unusual weather conditions (e.g., drought) (Estevez 2002; Palmer et al. 2015), hindered interpretation due to other aspects of naturally stochastic populations and interannual variability (e.g., favorable recruitment for specific taxa), and resulted in relatively speculative conclusions due to the lack of a suitable control.

Long-term monitoring programs in the Ten Thousand Islands area in southwest Florida, USA (Fig. 1) offer an excellent opportunity to evaluate watershed scale manipulation while addressing several of the issues outlined above. This area has multiple bays with contrasting freshwater alterations (i.e., enhanced, reduced, or relatively neutral compared with natural conditions) but otherwise similar physical properties. This region is part of the Greater Everglades Ecosystem and includes one of the largest expanses of mangrove forest in the USA. The maze of small mangrove islands, shallow bays, and passes here is home to large populations of wading birds and is a popular destination for recreational anglers seeking gamefish such as snook (Centropomus undecimalis), red drum (Sciaenops ocellatus), and mangrove snapper (Lutjanus griseus). Both the birds and gamefish are dependent upon the productivity of small prey fishes in the ecosystem (Bancroft et al. 1994; Patillo et al. 1997).

Fig. 1
figure 1

a Location of the Big Cypress Basin in SW Florida. b Locations of roads (gray) and canals (blue), and the trawl sampling grid and position of the water quality monitoring station (black drop symbol) in c Fakahatchee, d Faka Union, and e Pumpkin Bays. Average daily f salinity and g temperature among bays for 2000–2020 with error shading depicting the average minimum and maximum observed values on each day throughout the 20-year dataset

Full size image

These small-bodied and primarily juvenile prey-fish communities have been the subject of a long-term (~ 20 years) trawl-based monitoring program by the Rookery Bay National Estuarine Research Reserve (RBNERR). Paired with RBNERR’s water quality monitoring program, these long-term datasets provide an unprecedented opportunity to investigate the impacts of both higher and lower freshwater flow compared to a relatively natural system. The objectives in this study were to (1) quantify the differences in fish communities among seasons and bays with different hydrologies, (2) identify the primary species that are responsible for those differences based on their abundances, and (3) understand the relationships between structure of fish communities and their environment, specifically the influence of daily, weekly, and seasonal changes in salinity and temperature. Relationships identified from this study can shape expectations from watershed restoration actions in the region including the Comprehensive Everglades- and Picayune Strand Restoration Projects (U.S. Army Corps of Engineers 2004; Wingard and Lorenz 2014; South Florida Ecosystem Restoration Task Force 2020).

Methods

Study Area

The study area consisted of three estuarine bays in the Ten Thousand Islands region of southwest Florida: Fakahatchee, Pumpkin, and Faka Union (Fig. 1). These bays are similar in depth (~ 1 m on the sand/mud flats and ~ 5–7 m in channels during high tide), distance to the ocean (~ 6 km), tidal range (~ 1 m), and substrate (sand, mud, shell hash, and oyster bars), and all three are lined with a continuous fringe of red mangroves (Rhizophora mangle) (Carter et al. 1973; Colby et al. 1985; Browder et al. 1986; Shirley et al. 2004). Faka Union and Pumpkin are similar in size and have a single point source of tidal river input, whereas Fakahatchee is twice as large based on surface area and has 2 source rivers. The perimeter to area ratio is nearly identical for Faka Union and Pumpkin Bays (~ 0.006), but only half that value for the larger, more open Fakahatchee Bay (~ 0.003). Before human alteration of the watershed, these bays received freshwater from rainfall on the Big Cypress Basin that drained first as a slow sheetflow across the gently sloping landscape (~ 12 cm/km) and then gathered into seasonal streams and small tidal rivers that emptied into the bays (Carter et al. 1973). From the southwest, seawater from the Gulf of Mexico is brought into the bays on diurnal tides through the web of passes and channels in the Ten Thousand Islands (Fig. 1).

Over the last 100 years, however, watershed development and water control structures upstream from these bays have resulted in very different quantities of freshwater flow to the estuaries relative to historic conditions. Fakahatchee Bay retains the most natural flow regime, with sheet flow conveyed through the relatively undeveloped swamp called Fakahatchee Strand (~ 48,000 ha) (Yokel 1975; Shirley et al. 2004; Booth and Knight 2021) (Fig. 1b). This flow has been somewhat restricted by the construction of Highway US41 in the 1920s, logging roads in the 1940s, and region-wide depression of water tables, but water still passes southward through dozens of bridges and culverts, eventually flowing into Fakahatchee Bay through two small tidal rivers (Booth et al. 2014).

In contrast, the watershed for Faka Union Bay has been significantly altered through a failed suburban development (South Golden Gates Estates, SGGE) (Fig. 1b) in the 1960’s which included a massive network of 222 km of canals extending far up into the watershed. These canals effectively drained ~ 58,000 ha of the Big Cypress Basin through the Faka Union Canal (i.e., the channelized remains of the dredged Faka Union River) resulting in annual freshwater flow to Faka Union Bay ~ 10–100 times greater than in neighboring bays (Shirley et al. 2004; Booth and Soderqvist 2016). In fact, ~ 90% of the cumulative freshwater flow reaching the three study bays comes from the Faka Union Canal alone (Booth and Knight 2021). This was the case from the 1960s until 2021 when the canals were blocked and watershed restoration came online (Booth and Knight 2021). During the timeframe of this study, however, these canals quickly delivered large volumes of rainwater into Faka Union Bay instead of allowing its slower passage through natural sheet flow and a small tidal river.

In Pumpkin Bay, the opposite problem exists. Canals from SGGE, water control features, and agricultural development have resulted in much less freshwater reaching Pumpkin Bay which drains a smaller watershed (~ 5000 ha) than was present pre-development (Booth et al. 2014). In fact, only 1–2% of the total freshwater flow to the three study bays enters through Pumpkin River (Booth and Knight 2021). This has created a much more saline environment in the Pumpkin Bay system such that mean annual salinities are ~ 5 ppt higher in Pumpkin River than in Fakahatchee River, and ~ 10 ppt higher than in Faka Union Canal/River (Colby et al. 1985; Booth et al. 2014). High salinity conditions exceeding 35 ppt occur in Pumpkin River ~ 30–60% of the time depending on annual rainfall, whereas such conditions occur in Fakahatchee and Faka Union tributaries only ~ 10–30% of the year (Booth et al. 2014; Booth and Knight 2021). This difference in altered flow regimes among these three bays with otherwise similar physical environments provides an ecosystem-scale opportunity to investigate the effects of both more, and less, freshwater flow compared to natural conditions on community structure of estuarine fishes (Grange et al. 2000; Gillanders and Kingsford 2002).

The subtropical climate of southern Florida is characterized by significant differences in seasonal rainfall. Beginning in May or June of each year, cyclic changes in regional ocean temperature and atmospheric circulation result in almost daily afternoon showers and thunderstorms (Misra et al. 2018). This abrupt shift, which typically occurs over a matter of a few days, marks the onset of the early wet season (June–August) (Shirley et al. 2004) (Fig. 1). The end of this near-daily rainfall typically occurs in October (Misra et al. 2018) marking the late wet season (September–November) as drainage from accumulated rain in the watershed continues for some weeks after daily rains cease. Approximately two-thirds of the total annual rainfall occurs from June through October. Winter and spring months are marked by much less rainfall, declining watershed drainage, and increasing salinity in the estuaries. This period can be classified into the early-dry (December–February) and late-dry seasons (March–May). Maximum salinity, often in excess of 35 ppt, typically occurs in May in all three bays, and minimum salinity occurs in September (Christensen 1998).

The effects of water control structures in these watersheds are evident from seasonal salinity patterns (Fig. 1f). Pumpkin Bay remains more saline except during late dry season, when salinity peaks in all three bays. It experiences less dramatic declines in salinity during the wet season, with a typical minimum of ~ 20 ppt and exhibits an overall lower range in daily salinity relative to the other bays. Faka Union Bay experiences the greatest decline in salinity with the onset of the wet season, with a typical minimum of ~ 2 ppt. It experiences a greater range in daily salinities and remains fresher for longer into the dry season than the other bays. Fakahatchee Bay has a pattern intermediate between the other two bays.

Trawl Data

A long-term monitoring program targeting primarily small and juvenile fishes in this area was begun in 1998 by RBNERR (Shirley et al. 2004). Since 2000, four replicate trawls have been conducted monthly in each of the three study bays continuously for over 20 years with a pause from July 2013 to December 2015 due to funding constraints. Trawl sites were randomly selected using a 100-cell sampling frame overlaid on a nautical chart for each bay (Fig. 1c–e). Even though Fakahatchee is larger than the other bays, the size of the sampling frame is the same and it is positioned over the east side of the bay near the river outlets and away from the connection to Faka Union Bay where there is evidence that some canal flow enters Fakahatchee in the wet season (Browder et al. 1986). All sampling is conducted during a 4-day period each month and within 2 h of high tide. A 6 m wide otter trawl with 38 mm stretched mesh as the body of the trawl net and 3.2 mm bag liner was used for all sampling. Trawls were pulled into the current at a constant speed for ~ 0.19 km. All captured fish were identified to species or lowest possible taxonomic level, counted, measured (total length to nearest mm up to 20 individuals of each species), and then released. Contents of the four replicate trawls in each bay were combined into a single sample representative of the entire bay for each month. Counts for each species were summed and then divided by the total trawl distance (i.e., catch per unit effort [CPUE]) to standardize for variable trawl length and because some months did not have all four replicates completed (e.g., due to weather). Also, because presence of physical items on the bottom offer different structure and resources for fish and may influence their community structure, the bycatch volume of macroalgae, sponges, seagrass, mangrove propagules, and woody debris was also recorded.

Environmental Data

Estuarine species vary in their sensitivity to mean salinity and temperature as well as minima, maxima, and variation (Lorenz 1999; Faunce et al. 2004). In addition, fish communities can be influenced not only by the salinity and temperature observed at the time of the trawl, but also by the environmental conditions in the days and weeks leading up to each sample (Estevez 2002). To explore the most influential timeframe and relationship of these variables to fish communities, we calculated them at intervals of 1 day, 3 days, 1 week, 2 weeks, 1 month, 2 months, and 3 months before each trawl sample. Raw data for calculating these parameters were collected from automated water-quality monitoring stations that record data at 15-min intervals. These instruments are mounted on fixed poles in flow-through PVC housings near the primary freshwater input source for each bay (Fig. 1c–e) (for additional details see NOAA NERRS 2020). Mean and standard deviation were calculated from all water-quality readings during each time interval. Minimum and maximum are simply the very lowest and highest values recorded during each interval. Trawl samples with missing or erroneous environmental data due to sensor failure were excluded from consideration.

Statistical Analyses

Summary Statistics

We tested the effect of bay (Fakahatchee, Pumpkin, and Faka Union) and season (early wet, late wet, early dry or late dry) on total CPUE, species richness, and by-catch volume using an analysis of variance (ANOVA), followed by a Tukey multiple comparison test. Total CPUE was natural log-transformed and bycatch volume was natural log + 0.01 transformed to meet the assumptions of ANOVA. Only macroalgae and sponges were observed commonly enough for analysis of bycatch items. We also calculated mean CPUE for each species across each bay and season throughout the study period. All statistical analyses were conducted in R Version 3.6.1 (R Core Team 2019).

Fish Communities Compared Among Seasons Within Each Bay

First, we analyzed each bay separately, to understand how fish community structure differed among seasons (i.e., wet early, wet late, dry early, dry late) as related to changes in salinity and temperature. To prevent vagrants and extremely rare taxa from obscuring differences in fish communities among seasons, only species present in > 1% of the monthly sampling units in each bay were analyzed. Non-metric Multi-Dimensional Scaling (nMDS) was used on square-root transformed CPUE data using the metaMDS() function from the vegan package (Oksanen et al. 2019). More severe transformations and untransformed data were explored in preliminary analyses but yielded poor results such as unacceptably high stress in nMDS plots. In this analysis, a dissimilarity matrix comprised of all monthly samples within a bay was created based upon Bray–Curtis distances, and relationships among the samples were plotted in 2-dimensional space such that those with more similar fish communities are located closer together than samples with less similar communities.

To determine if fish communities differed significantly among seasons within each bay we used analysis of similarity (ANOSIM) with the anosim() function from the vegan package (Oksanen et al. 2019). If the overall ANOSIM was significant, then pairwise tests with a Bonferroni-adjusted p value were conducted to determine which seasons differed from each other. The resemblance statistic (R) was used to represent the size of the difference between communities (Clarke and Gorley 2006). When significant differences in community structure were detected among seasons, we used a similarity percentages analysis (SIMPER) (simper() in the vegan package, Oksanen et al. 2019) to determine which species were primarily responsible. This function calculates the proportion of community dissimilarity between seasons that is explained by each individual species.

Next, we sought to identify the environmental variables (i.e., minimum, maximum, mean, and standard deviation of salinity and temperature measured at various timescales from 1 day to 3 months prior to each trawl sample) that best related to the fish community in each sample. Specifically, we used the BIO-ENV function to find the set of environmental variables from each monthly sample whose dissimilarity matrix best correlated with the corresponding fish community dissimilarity matrix (Clarke and Ainsworth 1993; Oksanen et al. 2019). In this analysis, the relationship between the environmental and fish community matrices is expressed as a Spearman rank correlation coefficient. Due to the large number of correlated environmental variables, a two-stage process was used to determine which environmental variables best related to fish community structure. First, for each environmental parameter, we identified the single time interval (e.g., minimum salinity 1 day verses 1 week before a sample) that best correlated with the fish community. Then, considering only the best time-intervals, all possible combinations of a temperature and salinity variable (standardized) as well as individual salinity or temperature variables were tested. The environmental variable or variables with the highest Spearman rank correlation coefficient from the second stage of this analysis were deemed to best relate to the fish community.

Fish Communities Compared Among Bays Within Each Season

We then conducted a similar set of analyses to understand how fish communities compared among the bays during the four seasons. As before, species present in < 1% of the monthly sampling units (all bays) in each season were excluded, and nMDS was used on square-root transformed CPUE data. ANOSIM was used to determine if there was a significant difference in fish community structure among bays within each season, SIMPER was used to determine which species were primarily responsible for the differences among bays, and the two-stage BIO-ENV procedure was used to determine which salinity and temperature variables as well as time periods for measuring them (e.g., 1 day to 3 months prior to the sample) best correlated to community structure.

Results

Over the 20-year monitoring period (2000 to 2020), there were 1659 trawls with complete environmental data covering a combined total of 307 km. Pooling the monthly trawl samples in each bay resulted in 419 samples for analysis (146 in Fakahatchee, 165 in Faka Union, and 108 in Pumpkin) with differences in numbers among bays due to field logistics and maintenance of water quality stations. This large number of samples likely enabled high statistical power for detecting even small differences in fish communities among seasons and bays.

Over 260,000 individual fishes (mean TL of all fish = 52.4 mm ± 0.14 SE) were caught in the trawls representing 47 families, 76 genera, and 88 species or species groups (Electronic Supplementary Appendix A: length by individual species and CPUE by bay and season). Average total CPUE was significantly different among all three bays, with Pumpkin Bay having the highest, Fakahatchee having the lowest (30% lower than in Pumpkin), and Faka Union having intermediate values (Table 1). Mean species richness per monthly sample was significantly higher in Fakahatchee (~ 8 species) and Pumpkin (~ 8) than in Faka Union (~ 7). Both CPUE and species richness were significantly higher in the wet seasons than in the dry seasons. Mean CPUE was over twice as high in the wet seasons and richness increased from ~ 7 to ~ 9 species.

Table 1 Mean CPUE and species richness of fishes, mean volume of by-catch of algae and sponges, and number of monthly trawl samples for each bay and season (± SE)
Full size table

There was significantly less algae in Faka Union compared to either Fakahatchee or Pumpkin Bay, the latter two having ~ 4 times greater volume in the trawls (Table 1). Algal volume also varied significantly by season, with lowest volume in the late wet season and ~ 4 times higher values in the late dry season. There was significantly less volume of sponges in Faka Union than in either Pumpkin (~ 4 × higher) or Fakahatchee (~ 14 × higher than Faka Union). There were also significant differences in sponge volume by season, such that the late dry season had the smallest volume and the late wet season had ~ 3.5 times greater volume.

Fish Communities Compared Among Seasons Within Each Bay

For all three bays, nMDS and overall ANOSIM results indicated significant differences (p < 0.001) in fish community structure among seasons (Fig. 2). Stress values indicated acceptable ordinations. Dispersion was somewhat greater in the dry seasons. Pairwise tests for differences among seasons indicated that all seasons within bays were significantly different from each other (padj < 0.006); however, this is primarily the result of the large sample size and high power for detecting even small differences between groups. More informative are the R values, which convey the size of the differences between seasons. For interpretation, R values close to 1 indicate completely different fish communities between seasons whereas those close to zero indicate very similar communities (Clarke and Gorley 2006). In Fakahatchee, only the comparison between the late dry season and late wet season had a high R value (0.72), suggesting a large difference between fish communities. All other comparisons between seasons were below R = 0.4 indicating only moderate differences in community structure (Table 2), although the transitions between seasons showed a smooth cyclical pattern (i.e., early wet, late wet, early dry, late dry, and then back to early wet) in the nMDS (Fig. 2a). In Faka Union, the fish communities in the early and late wet season were quite different from those in the late dry season (R > 0.6; Table 2). There was also a moderate difference between the early wet and early dry seasons (R = 0.49). In Pumpkin Bay, the largest difference in fish communities was detected between the late dry and late wet seasons (R = 0.74; Table 2). All other pairwise comparisons between seasons showed only small to moderate differences (R = 0.13 to 0.42). Overall, these results demonstrate that the largest seasonal differences in fish community structure for all three bays occurred when comparing the late dry and late wet seasons.

Fig. 2
figure 2

nMDS plots comparing fish communities among seasons within a Fakahatchee, b Faka Union, and c Pumpkin Bay. Each point represents a pooled sample in each bay for each month in the 20-year dataset. Ellipses denote the standard deviation of each group centroid

Full size image
Table 2 ANOSIM R values from pairwise tests for differences in fish communities between seasons in Fakahatchee, Faka Union, and Pumpkin Bays. All comparisons were significantly different (padj < 0.006)
Full size table

The SIMPER analysis revealed which species were most responsible for the seasonal differences in the fish communities. Rather than showing all of the species in each comparison, only those cumulatively responsible for a majority (60%) of the differences in the seasonal fish communities are listed (Table 3a–c). The rest of the species contributed only small amounts (0–3%) to the differences between communities. In most seasons for all three bays, Mojarras and Anchovies were responsible for > 30% of the difference in fish communities. Pinfish (Lagodon rhomboides) also contributed to differences in communities between the early and late dry seasons in all three bays. The rest of the top species contributing to the differences in fish communities varied in their order of importance but always included blackcheek tonguefish (Symphurus plagiusa), silver perch (Bairdiella chrysoura), hardhead catfish (Ariopsis felis), lined sole (Achirus lineatus), and inshore lizardfish (Synodus foetens) for one or more seasonal comparisons in all three bays.

Table 3 Results of the SIMPER analysis showing the percentage of the difference between seasons explained by each species in (a) Fakahatchee, (b) Faka Union, and (c) Pumpkin Bay. Seasons (wet early, wet late, dry early, dry late) are abbreviated (e.g., dry late is DL)
Full size table

The seasons in Fakahatchee with the largest difference in fish community based on the ANOSIM were late dry versus late wet. Compared to late dry, the late wet season had 7 to 10 times more Mojarras (Eucinostomus spp.), Anchovies (Anchoa spp.), and Blackcheek Tonguefish, 16 times more Lane Snapper (Lutjanus synagris) and Lined Sole, and many more Catfishes (64 and 221 times more Hardhead and Gafftopsail [Bagre marinus]) (Fig. 3a). Not all species were more abundant. There were 10 times fewer Pinfish and 49 times fewer Pigfish (Orthopristis chrysoptera) in the late wet season compared to the late dry season.

Fig. 3
figure 3

Mean CPUE (log scale ± SE) by season for the species in a Fakahatchee, b Faka Union, and c Pumpkin Bay contributing to at least 60% of the difference in fish communities among seasons. Bars are left blank for those species that did not contribute to the 60% threshold within each bay

Full size image

In Faka Union, the seasons exhibiting large differences in fish community structure were late dry versus both the early and late wet seasons. The wet seasons had 4–17 times more mojarras, sand seatrout (Cynoscion arenarius), and clown (Microgobius gulosus) and green gobies (M. thalassinus), but 10–100 times fewer pinfish compared to the late dry season (Fig. 3b). In Pumpkin Bay, results were similar to those from Fakahatchee in that the late dry versus late wet seasons exhibited the largest difference in fish communities. In the late wet season, there were 5–51 times more mojarras, blackcheek tonguefish, lane snapper, hardhead catfish, lined sole, and sand seatrout compared to the late dry season, and 22 times fewer pinfish and no pigfish (Fig. 3c).

Spearman’s rank correlation (rs) values from the BIO-ENV procedure identified the time interval of salinity and temperature measurements that best related to the fish community structure in each bay (Table 4a–c). In Fakahatchee, mean salinity had its highest correlation with the fish community when measured over a time span of 8 weeks although all time periods were generally similar (Table 4a). Shorter timespans correlated best with the fish community for minimum, maximum, and standard deviation of salinity. Correlations for these were highest when measured at a time interval of 1 or 2 weeks, although shorter intervals of 1 to 3 days were correlated nearly as well. In contrast, all of the temperature variables (mean, minimum, maximum, and standard deviation) had highest correlations with the fish community in Fakahatchee when measured at longer time intervals of 8 to 12 weeks before the trawl samples. In Faka Union, compared to the other bays, the shortest time intervals of measuring salinity were the most related to fish communities (Table 4b). Mean and maximum salinity calculated over just the day prior to the trawl sample had highest correlation with the fish community. Minimum (1 week) and standard deviation of salinity (4 weeks) best related to fish communities in Faka Union at intermediate timescales. Similar to Fakahatchee, all temperature variables related best to fish communities in Faka Union when measured during the 8 to 12 weeks before the trawl samples. In Pumpkin Bay, only maximum salinity had highest correlation with the fish community when measured at a 1-week interval before the trawls. All other aspects of salinity and temperature had highest correlation to fish communities in Pumpkin Bay when measured at longer time spans of 2 to 3 months before trawl samples. Overall, this suggests that seasonal scales of temperature measurements best relate to fish communities in all three bays, whereas salinity is most related to fish communities at short time scales (days) in Faka Union, intermediate timescales (weeks) in Fakahatchee, and longer timescales (months) in Pumpkin.

Table 4 The Spearman rank correlation between the fish dissimilarity matrices for (a) Fakahatchee, (b) Faka Union, and (c) Pumpkin Bay, with salinity and temperature variables measured at different time intervals
Full size table

Considering only the best time intervals for each salinity and temperature variable using the highest Spearman rank correlation values from Table 4, we then tested each individually as well as all possible combinations of them to determine which dissimilarity matrices of variable (s) best related to the fish dissimilarity matrices for each bay. For both Fakahatchee and Faka Union Bays, the highest Spearman correlations (rs = 0.39 and 0.48, respectively) with the fish community were identified using only mean temperature measured over a time span of 12 weeks. No combinations of multiple temperature or salinity variables measured at any other time interval exhibited a higher correlation. In contrast, in Pumpkin Bay, the combination of variables with the highest correlation (rs = 0.43) to the fish community included minimum salinity and mean temperature recorded during the 12 weeks before the trawls.

Fish Communities Compared Among Bays Within Each Season

For all four seasons, nMDS and overall ANOSIM results indicated a significant difference (p < 0.001) in the community structure among bays, although overall R values were between 0.08 and 0.26, indicating only a small amount of actual dissimilarity among groups. Stress values indicated acceptable ordinations (Fig. 4).

Fig. 4
figure 4

nMDS plots comparing fish community structure among bays for the seasons a early wet, b late wet, c early dry, and d late dry. Each point represents a pooled sample in each bay for each month in the 20-year dataset. Ellipses denote the standard deviation of each group centroid

Full size image

Pairwise tests between bays revealed several significant differences in community structure among bays within each season; however, this may have been the result of the large sample size and high power for detecting even small differences between groups. The R values which report the actual size of the difference were more informative (Table 5). In the early wet season, the largest difference in fish communities was between Fakahatchee and Faka Union Bays. However, even this demonstrated only moderate differences in the fish community (R = 0.38). All other differences among bays in other seasons were either not significant or the actual differences were quite small (R = 0.1 to 0.25).

Table 5 R values from pairwise tests for differences in fish communities between bays in different seasons. All reported comparisons were significantly different (padj < 0.006)
Full size table

The SIMPER analysis identified which species are most responsible for the differences in the fish community structure among bays. Because only the early wet season showed even a moderate difference in community structure between Fakahatchee and Faka Union Bays based on the R value (Table 5), only SIMPER results for that comparison are presented. As was done for seasonal comparisons, only those species cumulatively responsible for a majority (60%) of the differences in the fish communities are reported (Table 6). Remaining species each contributed < 3% of the differences between communities. Mojarras were responsible for the largest share (18.8%) of the differences between communities, accounting for more than twice as much of the differences compared to the second most important species, anchovies (8% of the difference). Pinfish and silver perch each consistently accounted for 5–8% of the differences between bays.

Table 6 Results of the SIMPER analysis showing the percentage of the difference between Fakahatchee and Faka Union Bays that is explained by each species during the early wet season
Full size table

In the early wet season, Fakahatchee had 4–45 times more pinfish, silver perch, hardhead catfish, pigfish, code goby (Gobiosoma robustum), and gulf pipefish (Syngnathus scovelli) than Faka Union (Fig. 5). In contrast, Faka Union had 1.5–3.5 times more mojarra, anchovy, blackcheek tonguefish, and sand seatrout than Fakahatchee.

Fig. 5
figure 5

Mean CPUE (log scale ± SE) by bay for the species during the early wet season contributing to at least 60% of the difference in fish communities between bays based on SIMPER analysis

Full size image

Spearman rank correlation values from the BIO-ENV procedure including all bays combined indicated that salinity parameters were best related to fish community structure when measured at short time intervals (1–3 days before trawl samples) except for minimum salinity (4 weeks), although correlations for all time intervals were roughly equivalent for that variable (Table 7). All temperature variables best related to the fish community when measured at the longest time interval of 12 weeks.

Table 7 Spearman rank correlation between the fish dissimilarity matrix during the early wet season in all bays and salinity and temperature variables measured at different time intervals
Full size table

Considering only the time intervals for each salinity and temperature variable with the highest Spearman correlation values, we then tested all possible combinations of those variables to determine which exhibited the best relationship to the fish dissimilarity matrix during the early wet season. The combination of variables with the highest Spearman correlation (rs = 0.39) to the fish community included minimum salinity during the 4 weeks before the trawls and mean temperature during the 12 weeks before the trawls.

Discussion

As human populations in coastal watersheds continue to grow in the coming decades, understanding the impact of altered freshwater flow to estuaries will be a top information need for coastal managers who must balance diverging demands on freshwater resources (Doering et al. 2002; Gillanders and Kingsford 2002; Erwin 2009; Palmer et al. 2015). The Ten Thousand Islands provide a rare, ecosystem-scale opportunity to investigate the influence of both more, and less, freshwater flow compared to a bay with relatively natural hydrology. Here, we combined information from a 20-year trawl dataset and a water quality monitoring program to investigate how the difference in freshwater flow regimes between three bays with otherwise similar physical environments affects the community structure of estuarine fishes.

Despite major contrasts in salinity and watershed manipulation, the differences in fish community structure among study bays as well as seasons were not due to a different suite of species being present, but were instead the result of different abundances of a few common species. This occurred despite square-root transforming the abundance data in our analysis, which downplays the differences among highly abundant species and enhances the influence of less common ones. This finding is a robust confirmation of earlier, short-term assessments of the fish community in the area that found broad overlap in dominant species in these and other nearby estuaries along the Southwest Florida coast (Carter et al. 1973; Colby et al. 1985; Browder et al. 1986; Sheridan 1992; Shirley et al. 2004; Idelberger and Greenwood 2005; Tolley et al. 2006). Not only was a similar set of species responsible for seasonal differences in all three bays studied here, but more importantly, their relative abundance tracked with season in similar ways despite large differences in hydrology. This is one indication that season, rather than salinity, may be a stronger influence on community composition in this system.

In contrast, other estuaries with altered watersheds and seasonal fluctuations in salinity do experience greater change in species composition. For example, in some riverine estuaries in South Africa, emigration of some species during their wet season, and dramatic increases in abundance of anadromous fishes in the dry season, cause a more complete community turnover (Kanandjembo et al. 2001). In Puerto Rico, another riverine estuary that has been modified and loses virtually all of its freshwater flow during the dry season due to municipal diversion, experiences a major shift toward marine adapted species at the head of the estuary (Smith et al. 2008). Elsewhere in the Everglades ecosystem, in mangrove creeks with somewhat reverse hydrological conditions from the Ten Thousand Islands (i.e., low salinity in winter [< 5 ppt] and more saline [> 20–30 ppt] in summer due to the different manipulation of freshwater flow), several freshwater adapted species are found in abundance (Centrarchidae, Cyprinodontidae, and Cichlidae) during the low salinity winter (Faunce et al. 2004). The findings of the present study as well as those noted above, highlight the need to understand each system individually in the context of historical or natural conditions, seasonal cycles of migration and recruitment, and the severity and timing of changes in freshwater flow from water management (Alber 2002; Estevez 2002).

It was expected that there would have been at least a short-term seasonal shift toward freshwater adapted species in Faka Union due to canal flow; however, this was not detected. Even though Faka Union has salinity at least 50% lower than the other bays during the wet season, it still maintains an average salinity of 10 ppt. Faka Union also achieves high salinity conditions (36–40 ppt) annually during the late dry season as do the other two bays (Booth et al. 2014; Soderqvist and Patino 2010). This full range of salinity values every year, and lack of truly freshwater conditions at any time, may prevent the fish community of Faka Union from shifting to a more freshwater assemblage even briefly during the wet season.

Although differences in the fish community among bays were fewer than anticipated, some were still detected. For example, Fakahatchee and Pumpkin Bays had significantly different fish communities in the late dry season, when conditions in the bays are most similar, although the magnitude of the difference was very small. Only Fakahatchee and Faka Union had even moderately different fish communities, and only during the wet season when the abundances of many fishes are highest (Browder et al. 1986). The largest share of the difference between fish communities of Fakahatchee and Faka Union Bays was due to Mojarras, which accounted for more than twice as much of the difference compared to other species. Although Mojarras were common in both bays, they were by far more abundant in Faka Union. Sampling in 1982–1983 (~ 15 years after construction of Faka Union canal), Colby et al. (1985) found similar patterns between the two bays during the wet season for not only Mojarras, but also Pinfish, Silver Perch, and Pigfish, but the opposite pattern for Hardhead Catfish and some other species compared to the longer term dataset used here. Another difference of potential importance between bays was that there were significantly fewer species per trawl in Faka Union than the other two bays, which may translate into lower biodiversity overall (Browder et al. 1986).

Most of the species responsible for the differences in fish communities in this study are prey for gamefish such as spotted seatrout (C. nebulosus), red drum, snook, and tarpon (Megalops atlanticus) (Patillo et al. 1997). These prey species are tolerant of a wide range of salinity from near zero to hypersaline conditions, can be found over a diversity of habitats, and differ only in relative abundance among the bays (Yokel 1975; Browder et al. 1986; Patillo et al. 1997). The extent to which changes in the relative abundance of these species may affect economically-important gamefish via prey availability is unknown.

Spawning seasons for the dominant taxa and lag times for recruitment are also important to consider when interpreting the observed patterns. Timing of anchovy spawning may be more concentrated during fall months in South Florida waters (Powell et al. 1989). Silver perch and mojarras have a major spawning peak in spring and another possibly in the fall (Charles 1975; Patillo et al. 1997). For these species, recruitment of sizes susceptible to the trawl gear in the bays occurs a couple of months later. Pinfish and pigfish juveniles are known to arrive in estuaries during spring and summer (Yokel 1975; Patillo et al. 1997), and indeed both exhibited a large increase in our trawl samples in those months, followed by a large decline in early winter. This is also consistent with Browder et al. (1986) who observed more ichthyoplankton larvae during spring months, dominated by many of the same species that were highly abundant in the trawls. In a conclusion similar to other nearby estuaries (Idelberger and Greenwood 2005), it is these differences in spawning patterns that are likely the cause of the regular changes in fish communities among seasons, particularly the high abundance of small fishes caught during the wet seasons (Browder et al. 1986) following spring spawning.

What aspects of the environment best relate to fish communities? At first it, may seem that salinity could be the major driver since the largest difference in fish communities coincides with the largest seasonal difference in salinity (late dry vs. late wet) and salinity variation over the days and weeks before a trawl is related to community structure. However, longer term measurements of temperature (measured 2 to 3 months prior to the trawl samples) almost always had higher correlations with the fish community than salinity. More importantly, the final variable selection process indicated that short term patterns in salinity are relatively less important in this system compared to seasonal aspects of temperature. Even though these are highly modified watersheds with dramatically different salinity characteristics in the estuaries, a primary influence on fish communities appears to be simple seasonal patterns in spawning and recruitment of the dominant species which is best correlated with temperature measured over a scale of 2–3 months. This finding is consistent with observations in nearby natural systems. Seasonal differences in juvenile fish assemblages in two rivers elsewhere in Southwest Florida were due to spawning and recruitment periods rather than differences in environment (Idelberger and Greenwood 2005). Similar interpretations have been made in juvenile fishes in estuaries in Australia wherein salinity and temperature were not major contributors to differences in fish community among seasons and bays (Blaber and Blaber 1980). Instead, seasonal spawning patterns were suspected as the most important factors.

Results here are simply environmental correlations, not conclusive causes, and several studies have highlighted the challenges in disentangling the relative influence of correlated environmental variables. In South African estuaries, Kanandjembo et al. (2001) noted that in general, it is cooler and more saline in the winter and warmer and less saline in the summer which makes it difficult to distinguish which is a more important effect of the environment, or if patterns are due simply to the season when juveniles of these estuarine species recruit to the bays. Similarly, in Northeastern Florida Bay, Faunce et al. (2004) noted that the responses of many factors are correlated if not directly linked (i.e., salinity, depth, flora, fauna), preventing the identification of which variables are most responsible.

With correlations between the fish community, salinity, and temperature being relatively moderate in this and some similar studies (Blaber and Blaber 1980; Idelberger and Greenwood 2005), it is important to consider what other sources of variation may shape fish communities in these bays. For example, Fakahatchee Bay is a good control estuary for comparison due to its more natural freshwater flow, similar distance to the ocean, substrate types, and mangrove dominated perimeter; however, it is larger and slightly deeper, and has a smaller perimeter to area ratio than the other bays. These physical traits may contribute to differences in the fish community. However, note that in most cases, the fish community in Fakahatchee was more similar to Pumpkin Bay even though the landscape (e.g., bay size) of Pumpkin is more similar to Faka Union. Faka Union does have significantly less algae and sponges than the other two bays, probably due to the larger fluctuation of salinity in Faka Union (Butler et al. 2018). Oyster reefs also vary in health and living density among the bays (Volety et al. 2014; Loh et al. 2017). Such biological structures provide different habitats among bays (Bell 2008; van Lier et al. 2017) but are not suspected to be a major source of energy (McCarthy 2018). There are also different sources of carbon among the bays, with Faka Union receiving more in the form of detritus due to canal flow (Carter et al. 1973). Generally, however, the base of the food chain is particulate organic matter in the water column and sediments that are consumed by crustaceans, which are in turn consumed by small fishes, and ultimately by apex predators such as large gamefishes and sharks (McCarthy 2018). How much of the difference in fish community structure can be explained by each of these factors is difficult to determine because they are all correlated, largely constant throughout the sampling period, and cannot be experimentally manipulated at the scale of the bays. Instead, we may conclude that some combination of these variables may be responsible for a large share of the variability not explained by temperature and salinity.

The impact of altered flow on fishes in the Ten Thousand Islands has been intermittently investigated since the 1970s until the trawl monitoring program from this study began in 1998 (Carter et al. 1973; Yokel 1975; Colby et al. 1985; Browder et al. 1986; Shirley et al. 2004). Unfortunately, no data are available that quantified juvenile fish communities before, or in the immediate aftermath, of canal construction. These older assessments, while informative, had several inconsistencies that prevent direct comparisons to the dataset used here. These include differences in sampling gear, time of sampling (day vs. night, tidal stage, season), bays sampled, and sample size. Perhaps most importantly, previous studies were limited to a 1–3 year time series during which most were self-described as having environmental conditions that were atypical compared to long-term averages (e.g., a rainy dry-season). This is a problem for making rigorous inference in systems subject to large amounts of interannual environmental variability, as is the case in Southwest Florida (Idelberger and Greenwood 2005).

In comparison to these earlier works, the analysis here has several advantages for understanding environmental influences on community structure. Most importantly, this 20-year dataset encompasses multiple wet, dry, and normal years, as well as years affected by tropical storms and other rare events. This enables unprecedented long-term, ecosystem-level inferences that cannot be perceived by studies that last only a few years and are based upon a limited range of interannual variability. Also important, previous studies had to rely on time-of-sampling temperature and salinity values taken during trawls. In contrast, the data from continuous water quality monitoring stations that we used provides an understanding of the role of environmental variables over multiple timescales during the days and weeks before the trawls were taken (Faunce et al. 2004). Because of these advantages, we were able to demonstrate that longer term patterns of salinity and especially temperature are more influential than values observed on the day of the sample.

The hydrological problems of South Florida and the Everglades are not universal across the entire landscape. The effects of altered flow on fishes have been investigated extensively over several decades (Davis and Ogden [Eds.] 1994; Sklar and Browder 1998; Lorenz 2014) across a spectrum of habitats including uplands and headwaters (Kushlan 1980; Rehage and Loftus 2007), subtidal creeks (Faunce et al. 2004; Ley et al. 1994), canals (Serafy et al. 1997), channels and inlets (Browder et al.1986), and open bays (Montague and Ley 1993). Depending on proximity to drainage canals and position in the watershed, some locations receive too little fresh water and others too much or at the wrong time, compared to natural conditions. Location-specific recommendations are needed (e.g., Faunce et al. 2004) to appropriately restore hydrologic regimes in different parts of the landscape (Davis et al. 2005).

The state of Florida and US Federal Government have recognized the need to restore natural flow in the Big Cypress basin (South Florida Ecosystem Restoration Task Force 2020). Altered flow not only impacts estuarine fishes, but also hinders ecosystem services such as aquifer recharge and has negative effects on upland biodiversity, wading-bird foraging, and oyster reef and seagrass health (Carter et al. 1973; Bancroft et al. 1994; Davis et al. 2005; Volety et al. 2014; Loh et al. 2017). In response, government agencies implemented the Picayune Strand Restoration Project with the goal of “restor[ing] the historic hydroperiods and sheetflow patterns in the study area to the extent possible…” (U.S. Army Corps of Engineers 2004). To meet this objective, more culverts were installed under US Highway 41, roads from the defunct SGGE were bulldozed to fill drainage canals, and three massive pumping stations and water spreading basins were built to restore sheetflow southward across the watershed toward the Ten Thousand Islands (Booth et al. 2014; Booth and Soderqvist 2016). None of these individual steps have resulted in major changes in salinity at Fakahatchee, Faka Union, or Pumpkin Bay during the timeframe of this study (Booth and Soderqvist 2016). The linchpin of these integrated engineering components consists of plugging the Faka Union Canal north of Highway 41, which should enable these watershed restoration features to come fully online by 2024. Ongoing monitoring will be required to determine if the anticipated changes in hydrology to the study bays will be realized (Sklar and Browder 1998; Wingard and Lorenz 2014).

It should be acknowledged that based on the results here, major changes in the membership of the juvenile fish community in response to watershed restoration are not anticipated, although relative abundance may change. The range in flow characteristics of the three bays, and general similarity in fish communities suggests that conditions will remain within the tolerance of most fishes in all three bays and would continue to do so after watersheds are restored (i.e., less freshwater drained through Faka Union). Of note, however, despite expectations that Faka Union would have lower abundance based on previous assessments (Carter et al. 1973; Colby et al. 1985; Browder et al. 1986), the very robust dataset used here indicates that Fakahatchee actually has the lowest overall CPUE of all three study bays. Restoring and reducing flow to Faka Union to more closely match the hydrological seasonality of Fakahatchee may reduce the abundance of small fishes in the system overall, primarily due to fewer Mojarra and Anchovy, even though it may increase species richness in Faka Union due to increased availability of sponge and macroalgae habitat (Bell 2008; van Lier et al. 2017; Butler et al. 2018).

It is nevertheless important to continue monitoring these systems and revisit this trawl dataset in 3–5 years after the new hydrology has been in place to test this hypothesis and conclusively determine if the fish communities are changing with the restored flow patterns. Community response to any changes in managed flow will result in benefits to some species, detriment to others, and many will not have a response at all as long as the magnitude of the changes in environment are not far outside of the conditions experienced by the bays today. Provided that the annual flow includes both saline and freshwater periods, many of the species should continue to thrive. The amount of flow may just shift the relative abundance of the community one direction or the other in the estuaries either due to better survival or species shifting farther inshore or offshore to occupy preferred salinities. For example, in drowned river estuaries, greater freshwater flow can result in a repositioning of more saline-preferring community members down river during the wet seasons (Alber 2002; Smith et al. 2008; Palmer et al. 2015). However, moving seaward in the Ten Thousand Islands means moving into a different habitat. Instead of the shallow bays that were the focus of this study, moving seaward means a landscape dominated by narrow channels and a much higher ratio of mangrove fringe to open water. Faster currents and more mangrove-edge associated predators could reshape relative abundance of fishes in less predictable ways.

This study focused on relative abundance of the fish communities in the bays as a whole, but there are many other research opportunities possible from this large and long term dataset. Individual species with more narrowly defined environmental tolerance may be better indicators of environmental change (Alber 2002; Doering et al. 2002; Estevez 2002; Palmer et al. 2015). Fish biomass, spatial distribution within bays (e.g., channel, flats, fringe), and changes in the time series such as gradual shifts due to climate change (Ross et al. 2000; Erwin 2009; Krauss et al. 2011; Michot et al. 2017) or in response to discrete events such as hurricanes (Radabaugh et al. 2020) or drought years, also represent topics that may offer insight into the relationships between fishes and watershed condition, stress, or restoration. Indeed, it will be important to continue monitoring both the estuarine fishes as well as abiotic variables as climate impacts continue and potentially shift local environmental conditions beyond those experienced during the present time span of this dataset.