Abstract Climate change is leading to significant alterations to ecosystems all over the world and some of the resulting impacts on fish and fisheries are now becoming apparent. Estuaries, which are highly susceptible to climate change because they are relatively shallow and in close proximity to anthropogenic stressors, provide habitat to many fish species at a critical time in the life history, after transport and just prior to settlement in nurseries. Despite this, the long-term impacts of climate change on larval fish at this critical location/stage in the life history are not well documented. The larval fish assemblage of a coastal estuary was sampled once per week for twenty-six years at a fixed location in southern New Jersey, USA. We used ordination and regression analysis to evaluate the whole assemblage, individual species/family occurrence, and trends in total density and diversity over that time. The larval fish assemblage changed significantly in response to warming water temperatures. In addition, approximately one quarter of the species/families in the assemblage exhibited a statistically significant trend in individual occurrence over time. Of these, all five of the five northern-affiliated species decreased in occurrence while 18 of 21 southern-affiliated species increased in occurrence. Finally, total fish density and species diversity increased over the course of the study. The non-uniform response of the species/families in this larval assemblage is similar to what has been documented in other studies that evaluated the temporal trend of open ocean juvenile and adult fish assemblages.

Citation: Morson JM, Grothues T, Able KW (2019) Change in larval fish assemblage in a USA east coast estuary estimated from twenty-six years of fixed weekly sampling. PLoS ONE 14(10): e0224157. https://doi.org/10.1371/journal.pone.0224157 Editor: Heather M. Patterson, Australian Bureau of Agricultural and Resource Economics and Sciences, AUSTRALIA Received: June 27, 2019; Accepted: October 7, 2019; Published: October 23, 2019 Copyright: © 2019 Morson et al. This is an open access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited. Data Availability: The data underlying the results presented in the study are available from: https://www2.dnr.sc.gov/seamap/. Funding: The collection and analysis of the data for this long time-series came from multiple and disparate sources over the 26 years of this study. The major internal sources at Rutgers University included the Center for Coastal and Environmental Studies, the Institute of Marine and Coastal Sciences and the Department of Marine and Coastal Sciences. External sources included NOAA Sea Grant, the Electric Power Research Institute, the Jacques Cousteau National Estuarine Research Reserve and New Jersey Department of Environmental Protection to KWA and TG. The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript. Competing interests: The authors have declared that no competing interests exist.

Introduction Climate change is a topic of increasing importance because of the current and future implications of its effects, especially in marine and estuarine ecosystems [1– 4]. Increased temperatures, shifting winds, rising water levels, intensified storms and changes in pH and ocean currents [5–7] are only some of the effects of a changing ocean. These effects are leading to significant alterations in coastal ecosystems. For example, numerous works now demonstrate how climate change has impacted fish populations [8–15] and associated fisheries [16–21]. Estuaries, which are important as nurseries for both larval and juvenile stage fish [22–24], are highly susceptible to climate change because they are relatively shallow and thus have low thermal inertia, and generally are in close proximity to anthropogenic stressors [25–28]. Evidence of increasing temperatures and large fluctuations in salinity due to climate change has been documented in estuarine systems in the northern and southern hemispheres [25, 29–32]. This includes along the east coast of the U.S. in large estuaries, like the Hudson River estuary and Narragansett and Delaware Bays, and small ones, such as the Mullica River-Great Bay estuary [23, 33, 34]. Fish larvae that are transported to estuaries may therefore be even more affected by climate change than larger juveniles and adults that have the ability to move to more favorable areas. Long-term ecological monitoring programs are critical to resolving the issues associated with climate change [e.g., 8, 35–38]. The Rutgers University Marine Field Station (RUMFS) in southern New Jersey has been conducting larval fish sampling weekly at one place in a coastal estuary, behind Little Egg Inlet, under a fixed protocol for 26 years. To date, the assemblage at this site appears similar to adjacent inlets and estuarine thoroughfares and thus is representative of the region [39]. The larval fish assemblage in this estuary was previously described using the first six years of this time series [40]. In addition, it has been compared with samples collected from other estuaries to gauge regional synchrony or to evaluate life history for specific species such as Anguilla rostrata [41, 42], Conger oceanicus [43], Paralichthys dentatus [44, 45], Brevoortia tyrannus [46], Micropogonias undulatus [47], and Pseudopleuronectes americanus [48, 49]. Synthesis of research to date suggests that changes in abundance of selected species are the result of climate change, recovery of local spawning stocks, or other factors that could only be detected with a long time series [23, 47, 49], but no work has yet evaluated the longer-term temporal trends in the entire larval fish assemblage for this data set. Our specific objective for this paper is to document change in the assemblage of larval fishes in this estuary as estimated from twenty-six years of weekly sampling. Evaluating change at this critical stage in the life cycle and at this critical estuarine location in the coastal ocean system is important for two reasons. First, measuring change at the larval stage addresses distribution at an important, well-defined point in the reproductive cycle because larval abundance is affected by events of the very recent past. The duration of competency for larval settlement is short, typically days. Thus, response signals (fluctuations in larval abundance over time) are less likely to be obfuscated by effects of mixed age, maturity, diet, habitat choice, and other ontological or size correlates expected of a sample of multiple year classes as it would be in fishery-dependent trawl samples, though there is a possibility that these effects could be lagged and thus still manifest in the larval signals. Second, interannual variation in the assemblage of estuarine larvae of open ocean spawners are likely to be more highly correlated with year-class strength than similar metrics calculated from open ocean larval surveys. This is because a large part of the critical phase from spawning to larval settlement, including coastal transport, is complete when larvae arrive in coastal estuaries. This accounts for the influence of, for example, match-mismatch of larvae with predators or food resources at those earlier stages.

Methods Study site Larval fish sampling was conducted behind Little Egg Inlet, New Jersey in the northeast United States continental shelf ecosystem (Fig 1). Little Egg Inlet, a relatively unaltered inlet, is the primary source of Atlantic Ocean water that enters two estuaries, a drowned river valley (Mullica River-Great Bay) and an adjacent barrier beach lagoon (Barnegat Bay), and is split by the Sheepshead Meadows peninsula. This dual estuary has a broad, seasonal temperature range (-1.2° to 35° C; [50]), a moderate tidal range, and an average depth of 1.7 m. Sampling was conducted from a bridge over Little Sheepshead Creek (water depth ~3m), a thoroughfare connecting Great Bay and Barnegat Bay across Sheepshead Meadows located 2.5 km from Little Egg Inlet (Fig 1). Atlantic Ocean water flows into the estuary through Little Egg Inlet during flood tides, and portions are diverted into the mouth of Little Sheepshead Creek [48, 51]. Larval fish samples collected from this location are representative of dynamics occurring in the estuary proper [e.g., 31, 40, 48], at other local inlets and thoroughfares [39], as well as the adjacent inner continental shelf [45, 52]. PPT PowerPoint slide

PowerPoint slide PNG larger image

larger image TIFF original image Download: Fig 1. Larval fish sampling location behind Little Egg Inlet in New Jersey (inset). https://doi.org/10.1371/journal.pone.0224157.g001 Sampling methodology Larval fish entering the estuary were collected with a 1-m diameter, circular plankton net (1-mm mesh). Three 30-minute deployments at 1.5 m water depth during night-time (to reduce gear avoidance) flood tides were made weekly from 1990 to 2015. A General Oceanics flowmeter was attached to the mouth of the net to gauge flow for calculating the water volume sampled for each deployment. Surface water temperature and salinity were recorded at the beginning and end of each sample set. Larval fish were preserved in 95% ethanol and subsequently identified to family or species and counted. Up to 20 individuals from each species for each deployment were selected at random and measured and staged as preflexion, postflexion, or flexion. Analytical approach Larval fish catch data were organized and analyzed at two levels, species and family, and all analyses were repeated and reported for both levels. In addition, the assemblage was evaluated for the subset of data that included those species or families that occurred in at least 1% of the samples over the time series, hereby referred to as commonly occurring, while analyses of individual species or families were performed on the entire data set. Running analyses at both the family and species levels allowed us to make comparisons at these levels as has previously been done for these taxonomically difficult larvae [59]. The mean density (fish/1000 m3) of each species or family for each weekly sampling event was calculated from the three replicate tows. The annual mean density for each species was then calculated by averaging the date-specific means. This approach resulted in a year-by-species matrix of annual mean density for each species or family. To identify significant changes in species composition over time, we used this matrix to evaluate change in total larval density and species diversity, change in assemblage for commonly occurring species and families, and change in occurrence for all individual species or families. This multifaceted approach helped to identify recurring patterns in the data at the individual and assemblage levels that are robust to different analytical methods and are therefore likely to indicate true change in the larval fish assemblage. Analyses applied R v3.3.1 statistical computing software [53]. Whole assemblage analysis. For commonly occurring species and families, the catch matrices described above were evaluated with principle components analysis (PCA) using the dudi.pca function in the ade4 R package [54]. To help equally weight the contribution of low catch and high catch species, data were centered and standardized before the PCA analysis. A statistical approach was used to determine dimensionality for the PCA to reduce the chance of interpreting noise in the data or of missing significant trends [55]. To estimate the total number of principle components to retain for interpretation, we used the testdim function in the ade4 R package [54]. This function calculates the RV coefficient for matrices with alternative dimensions to evaluate whether additional dimensions add relevant information [55]. For each principle component retained, we then regressed the component scores across year and used the lm function in the base R package to evaluate whether there was a statistically significant trend over time. Individual/Family-level analysis. Some species or families may occur at any given time in a sample at the tail ends of their environmental niche, and therefore may only be present intermittently when conditions are most favorable. As described previously, these rare (observed in <1% of samples) species or families were removed from the whole assemblage analyses. However, a trend in individual occurrence in samples from a high temporal-resolution, fixed sampling location could indicate a shift in the geographic range minima or maxima for a given species or family. Therefore, we also evaluated the change in presence/absence over time for all species and families. We used a generalized linear model with a logit-link function and a binomial error distribution to evaluate the change in probability of capture over time for each species. Logistic regression models were fit using the standard glm function in the base R package [53]. We assigned a binary response (presence/absence) to each species for each sampling event and regressed this response across year. Those species or families that exhibited significant temporal change (p<0.05) were then categorized as declining or increasing in occurrence over time depending on whether the slope of the best fit line was negative or positive, respectively. In addition, some species and families were categorized as either northern (those typically spawning north of Cape Cod, Massachusetts, including Georges Bank) or southern (those typically spawning south of Cape Hatteras, North Carolina) and local (those typically spawning in the Middle Atlantic Bight) based on Able and Fahay [23] and references therein. For this subset of species, we identified how many exhibited statistically significant increasing or declining trends in occurrence. Total density and diversity analysis. To evaluate whether there was a net change in total density or species diversity over time, we first calculated the Shannon and Simpson indices of diversity using the diversity function in the vegan R package [56]. Total density and each of the diversity indices were then independently regressed across year using the lm function in the base R package to evaluate whether there was a statistically significant trend over time [53].

Discussion Long-term monitoring programs are critical to our understanding of how ecosystems are responding to climate change and climate related impacts. Programs that document larval fish assemblages through time, in particular, could serve as early indicators of ecosystem change or of shifting distributions because larval fish occur at low trophic levels [57–59]. Furthermore, documenting change in estuarine larvae from open ocean spawners at the critical phase in development just prior to settlement may provide a more reliable index compared to, for example, open ocean larval fish where larvae have not yet been subjected to coastal transport and associated match-mismatch with resources, predators and prey. However, much of the work documenting the impact of climate change on fish community composition has focused on adult assemblages [e.g. 10, 12, 14]. When larval fish assemblages have been evaluated, the time series have either been short in duration [e.g. 60, 61, 62], focused on larval assemblages from the open ocean [e.g. 4, 58, 63], or when focused on estuarine larvae, evaluated a specific species [e.g. 41, 45, 47] or a selection of species [22, 31]. To our knowledge this is the first time a long-term time series of an entire larval fish assemblage in a coastal estuary has been evaluated for a climate-related trend. Enhanced warming of the Northwest Atlantic Ocean and Shelf is anticipated to have major consequences for the marine ecosystems there. With this work we identified four significant temporal trends that provide important insight for how larval fish assemblages, for both transient and resident species, in coastal estuaries are responding to rapid environmental change. First, the overall composition of larval fish at the sampling site changed over time and this was true whether the assemblage was evaluated at the finer, species level or at the coarser, family level. Therefore, this adds to a growing body of literature demonstrating that larval, juvenile, and adult fish communities are responding rapidly to changes in the environment across the world [e.g. 12, 58], including in the Northwest Atlantic Ocean [4, 10, 14, 31]. This approach also helps to address the question of taxonomic sufficiency for difficult larval fishes [59]. We did not attempt to identify the mechanisms controlling the observed change in larval assemblage, partly because the larval fish in this assemblage come from wide ranging origins in the ocean and estuary and are therefore influenced by a range of local and regional conditions that were not measured. However, if we consider the first two PCA axes as climate trends, it suggests 25% of the change in the larval assemblage could be attributed to changing climate. Previous work has demonstrated that bottom temperature can explain more than half of the variation in adult assemblage shifts for species collected in open ocean trawl surveys [12]. Second, when we evaluated change in occurrence at the individual species and family levels we found that only approximately one quarter of the species or families in the larval fish assemblage exhibited statistically significant trends in occurrence and that the significant trends were not unidirectional. Previous evidence has demonstrated that species within an assemblage can have varying responses to changes in the environment [4, 12, 64, 65] and recent work has demonstrated that these responses by fish, as poikilotherms, are likely influenced by a thermal preference specific to each species [66]. Given the diversity of taxa found in the larval fish assemblage evaluated with this work, it is not surprising that some species exhibit no change in occurrence, while others are increasing or decreasing in occurrence. If changes in larval fish assemblage are reliable harbingers of change in juvenile and adult assemblages, these non-uniform responses of individual species within the larval community are likely to alter trophic interactions of adults [67, 68]. However, the overall ecosystem impacts of these changes may be buffered when interacting species are replaced by those with a more tolerant thermal range that also occupy similar trophic guilds [69]. For instance, we observed a decline in the occurrence of Menticirrhus spp. and a concurrent increase in the occurrence of M. undulatus, species at similar trophic levels. Given the changes in the larval fish assemblage we observed, future efforts could focus on modeling how these changes are likely to impact interactions between species in this assemblage or in their subsequent recruits and what those changing interactions suggest about supply to nursery habitats. Third, for nearly every species or family that exhibited a significant trend in occurrence that was also categorized as either southern or northern relative to the fixed sampling location, the trend was as expected given warming conditions [70] and given previous evidence that larvae and adult fish in this region are shifting predominately northward [4, 12]. That is, northern species, or those likely to have a thermal range skewed toward colder water temperatures, declined in occurrence, and southern species, or those likely to have a thermal range skewed toward warmer water, increased in occurrence, over the same time frame that the region experienced rapid warming. A similar trend was found when a shorter section of this time series was evaluated for a subset of species in the assemblage [31]. Fourth, over 70% of the species that exhibited significant change over the time series increased in occurrence, resulting in an increase in total diversity, and likely total abundance. This is most obvious for the larvae of species that occurred for the first time late in the time series and were consistently present thereafter. For example, C. bosquianus was not observed in the assemblage from 1990 to 2006, but was observed for six of the following nine years. Similarly, G. oceanicus was not observed from 1990 to 2004, but was observed every year after (2005–2015). As a result, these species have likely become residents of this estuary over the duration (twenty-six years) of the times series we evaluated. A trend of increasing diversity is supported by previous work [71], including for juvenile and adult fishes in the same region [12, 66], and suggests either a range expansion or a range shift is occurring for some species. When trends in marine species richness were evaluated for trawl data collected from nine different regions around the world (Gulf of Mexico, Northeast US, Eastern Bering Sea, Southeast US, Gulf of Alaska, West Coast US, Newfoundland, Aleutian Islands, and Scotian Shelf), eight showed positive trends and four had statistically significant trends [72]. The authors offer that an expansion in the range of transient species over time could be causing a net increase in richness. The impact of shifts in fish distribution, like the ones documented here, extends beyond the ecological consequences for the larval community and the subsequent recruits and adult fish. Shifts in marine species distributions also present a challenging set of obstacles for how fisheries are managed because management has traditionally relied on political boundaries for surveying populations, assessing stocks, and allocating quota [20, 21, 73]. These obstacles will require unique and innovative solutions that take into account alternative management scenarios under different climate projections [74, 75]. Time series of larval fish assemblages like this work could prove useful in this endeavor because observing and understanding the delivery of larval fishes to estuaries, both in terms of abundance and timing, provides an opportunity for predicting change [76]. Efforts are already underway to predict changes in distribution of marine species using habitat projection models [e.g. 66, 77], however these models can be highly uncertain [70]. Perhaps data collected from larval fish time series could be used in tandem with habitat projection models to better predict shifts in distributions of juveniles and adults before they happen, and this knowledge could serve as a tool for fishery managers to develop proactive solutions to management problems before they occur.

Acknowledgments We would like to thank numerous technicians at Rutgers University Marine Field Station, especially T. Malatesta, R. Hagan, J. Rackovan, J. Caridad, M. Shaw, R. Larum, S. VanMorter, and C. Denisevich, and the Jacques Cousteau National Estuarine Research Reserve volunteers, especially P. Filardi, S. Zeck, T. Siciliano, T. Bonovolanta, and E. Lesher for assistance with data collection and entry. R. Hagan organized the synoptic collection and sorting of larvae during most of these efforts. Work was approved by Rutgers University under IACUC Animal Protocol 88–042.