Change in larval fish assemblage in a USA east coast estuary estimated from twenty-six years of fixed weekly sampling

Summary

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.

Keywords:

Species diversity – Oceans – Larvae – Marine fish – Estuaries – Eels – Inlets – Anthropogenic climate change

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].

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 m³) 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].

Results

Sample description

A total of 1,252 weekly samples were analyzed over the course of twenty-six years from 1990 to 2015. Over 650,000 fish larvae were collected in these samples and identified to family or species. Representative fish from sixty families were observed and forty were observed in at least 1% of the samples (Table 1). A total of 138 different species were identified with sixty-five of those present in at least 1% of the samples (Table 1).

Whole assemblage response

Following PCA of the commonly-occurring species, a single component was retained for subsequent analysis, while for the commonly-occurring families analysis, two axes were retained (Table 2). The single axis retained for analysis of the species data accounted for 14% of the variation in the catch, while the two axes retained for analysis of the family data accounted for 14% and 12% of the variation in the catch, respectively (Table 2).

There was a significant (p<0.0001) annual trend in the first principle component for the species-level PCA (Fig 2A). Some species increased their relative contribution to the assemblage over time. These species, all of which loaded heavily on the first component in the positive direction, included Engraulidae sp., Bairdella chrysoura, Clupeiformes sp., Gobiosoma sp., Microgobius thalassinus, Clupeidae sp., Lagodon rhomboides, Sciaenidae sp., Gobiosoma ginsburgi, and Leiostomus xanthurus. Other species decreased their relative contribution to the assemblage over the course of the time series. These species all loaded heavily on the first component in the negative direction, and included Anguilla rostrata, Clupea harengus, Enchelyopus cimbrius, Gasterosteus aculeatus, and Ammodytes sp.

For the family-level PCA, both the first (p = 0.0005) and second (p = 0.0417) principle components exhibited significant annual trends (Fig 2B and 2C). For species weighting heavily on the first principle component, paralichthyids, clupeids, ophichthids, batrachoidids, serranids, elopids, and atherinopsids all increased their relative contribution to the assemblage over the time series, while ammodytids, gasterosteids, lotids, ophidiids, fundulids, anguilids, pholids, and belonids decreased their relative contribution (Fig 2B). Of species weighting heavily on the second principle component, engraulids, muglids, eleotrids, and gobiids increased their relative contribution to the assemblage over the time series, while syngnathids, scophthalmids, atherinopsids, uranoscopids, anguillids, and scombrids decreased their relative contribution (Fig 2C).

Individual/Family-level response

Thirty-six (26%) out of 138 species and sixteen (27%) out of sixty families exhibited a statistically significant trend in occurrence over the time series (Table 3, Figs 3 and 4). Of those, ten families increased in occurrence while six declined; twenty-eight species increased in occurrence while only eight declined (Table 3, Figs 3 and 4). While only approximately one quarter of the entire assemblage changed in occurrence, for the commonly-occurring species in the assemblage (those observed in >1% of the samples), twenty-nine (45%) out of sixty-five exhibited a statistically significant trend in occurrence (Table 1, Fig 3).

Some species that showed statistically significant declines in occurrence over time were also identified in the PCA analysis as contributing relatively less over time to the species composition (Figs 2A and 3). These included Clupea harengus, Gasterosteus aculeatus, and Ammodytes sp. Similarly, most species that showed patterns of increasing occurrence were also identified in the PCA analysis as contributing more to the species composition over time, including Engraulidae sp., B. chrysoura, Clupeiformes sp., Gobiosoma sp., M. thalassinus, Clupeidae sp., L. rhomboides, Sciaenidae sp., and G. ginsburgi (Figs 2A and 3).

A similar trend was observed for families that showed significant changes in occurrence. Paralichthyids, clupeids, serranids, elopids, eleotrids, and gobiids both increased in occurrence and increased in their relative contribution to the assemblage, while ammodytids, ophidids, fundulids, and atherinopsids both decreased in occurrence and decreased in their relative contribution to the assemblage (Figs 2B, 2C and 4).

For those species/families that showed significant change in occurrence over time that were also categorized as northern species, five out of five (100%) declined in occurrence, including ammodytids, Ammodytes sp., C. harengus, G. aculeatus, and stichaeids (Table 3; Figs 3 and 4). Out of twenty-one species/families that showed significant change in occurrence that were also categorized as southern species, eighteen (86%) exhibited increasing occurrence. These included B. chrysoura, Chasmodes bosquianus, Ctenogobius boleosoma, Dormitator maculatus, eleotrids, elopids, Elops saurus, Engraulidae sp., Erotelis smaragdus, Gobionellus oceanicus, Hyporhamphus meeki, L. rhomboides, Megalops atlanticus, M. thalassinus, M. undulatus, Opisthonema oglinum, and Symphurus plagiusa (Table 3; Figs 3 and 4).

Total density and diversity response

Species diversity and total fish density increased over time (Figs 5 and 6). For the Shannon index there was a statistically significant (p = 0.0244) trend at the alpha = 0.05 level, while the Simpson index and total density trends were only significant (p = 0.0998 and p = 0.0783, respectively) at the alpha = 0.10 level.

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.