Relationship between fish assemblage structure and predictors related to estuarine productivity in shallow habitats of a Neotropical estuary

The variability of fish assemblage structure with respect to seasonality in salinity and productivity remains to be elucidate to many Neotropical estuaries. In this study, we hypothesized that salinity gradient and a set of variables related to ecosystem productivity drive community parameters in the shallow-water fish assemblage of the north-south axis of the Paranaguá Estuarine Complex (Southern Brazilian coast). Samples were taken with beach seine monthly from May 2000 to April 2001. Supporting our hypothesis, richness and abundance increased with turbidity, warmer waters of the rainier summer seasons, which are more productive. This environmental setting favors reproduction, as well as juvenile recruitment and growth, whose intensities are highest in this period. Highest abundance was found in inner areas, which may be explained by greater food and habitat availability. Richness was higher in more saline waters, due to the proximity of the rich pool of marine fish species. We suggest that local human interventions ( e.g. , dredging) should be avoided during the rainy seasons that are critical for species life cycles. Salinization, low estuarine productivity, and warmer waters, which are expected with climate change and human impacts in the local watershed, could affect the integrity of the local fish assemblage.


INTRODUCTION
Global climate change has significantly impacted marine and estuarine ecosystems around the world (Harley et al., 2006). These impacts induce changes in local climatology, such as temperature, freshwater flow, winds, currents, and tide-induced circulation, in addition to changes in water quality and biodiversity of estuarine systems (Garcia et al., 2003;Martinho et al., 2009). There is currently a global concern regarding the loss of biodiversity, especially when related to anthropogenic pressures that impact ecosystem goods and services (Cardinale et al., 2012).
Estuaries are physical and biological-transition environments between the continent, rivers, and the ocean. They have significant ecological importance and have been pointed as one of the most productive ecosystems on the planet (Costanza et al., 1997). Estuaries are also highly relevant to other contiguous environments, serving as, for example, nurseries for several marine species and as sources of nutrients to shelf habitats. They are naturally susceptible to large environmental fluctuations, which lead to changes in the communities of primary producers and, consequently, in fish fauna (Blaber, 2000;Mclusky, Elliott, 2004;Able, 2005).
Fishes dominate the estuarine nekton and have been widely used as indicators of environmental changes. Given the dynamics of estuaries, fish communities are generally characterized by relatively low species richness that are well adapted to such conditions, and are limited by the range of strategies used to deal with environmental fluctuations. Thus, these assemblages usually vary widely in space and time (McLusky, Elliott, 2007;Martinho et al., 2009).
3/20 ni.bio.br | scielo.br/ni Luís Henrique M. C. Vergès, Riguel F. Contente, Camila Marion, Cívil P. C. del Castillo, Henry L. Spach et al. Shallow estuary areas such as tidal flats and estuarine beaches are vital for fish early development. In these environments, they find abundant food and protection from predators, which are more abundant and actives in the deepest regions of the estuary. Consequently, an important part of the fish assemblages in estuarine shallow areas is composed by juveniles (Blaber, 2000;Elliott, Hemingway, 2002;Mclusky, Elliott, 2004;Able, 2005).
Several factors are related to the structure of shallow-water fish fauna in estuaries. Biological interactions such as competition and predation may be involved in such assemblage structuring processes (Elliott, Hemingway, 2002), as well as abiotic factors, such as temperature, salinity, water transparency (Aguirre-León et al., 2014), dissolved oxygen , climatic conditions (Garcia et al., 2003), and oceanographic conditions (Miranda et al., 2002). Salinity and temperature have been recorded as the main abiotic factors regulating fish communities in terms of species richness and abundance, respectively (Whitfield, Elliott, 2002).
Many aspects of the role of abiotic factors on the richness and abundance of shallow water fish faunas in Neotropical estuaries still need to be better elucidated. A remarkable is the behavior of richness and abundance with respect to the conspicuous seasonal variability of the estuarine hydrological conditions (Barletta et al., 2003;Blaber, Barletta, 2016). In the rainy periods, a strong continental input of freshwater and nutrients reach the estuary. This, coupled with higher temperatures, increases substantially the local primary and secondary production rates. Such favorable conditions increase fish influx from the ocean to the estuary, increasing estuarine fish richness and abundance. Many species reproduce (or increase their reproductive activity) during conditions of high food supply to maximize juvenile survival. In dry periods, an opposite pattern occurs when lower continental input and lower temperature imply in lower productivity and, consequently, lower rates of fish colonization (i.e., lower richness and abundance) in the estuary (Blaber, 2000;Elliott, Hemingway, 2002;Barletta et al., 2003;Mclusky, Elliott, 2004;Barletta et al., 2005). Such seasonal model has not yet been objectively tested in many Neotropical estuaries, such as the Paranaguá Bay Estuarine Complex -PEC, southern cost of Brazil.
Statistical models that correlate fish responses to environmental conditions have been widely used to investigate the roles of estuarine environmental factors and hydrology on fish distribution, occurrence, and abundance (Nicolas et al., 2010;França et al., 2011). Currently, the use of models to predict structure and functioning of ecosystems has been shown extremely efficient (França, Cabral, 2015). In studies involving coastal and marine ecosystems, and with regard specifically to modeling of fish richness and abundance, generalized linear models (GLMs) has been used more predominantly (Nicolas et al., 2010;França et al., 2011;França, Cabral, 2015). França et al. (2011) used GLMs to identify the role of temperature and salinity as the most important for variation in fish species richness along the Portuguese coast. Gaussian linear models, when their assumptions are met, may also be appropriated tools to model biodiversitypredictors relationships (Ives, 2015). In this study, we used statistical models to test the hypothesis that a set of variables related to biological productivity and salinity gradient of the estuary drive the total abundance, richness, and structure (composition and abundance of species) in the shallow-water fish assemblage of the north-south axis of the subtropical Estuarine Complex of Paranaguá Bay (Southern coast of Brazil).

Study area.
The study area is located in the Paranaguá Estuarine Complex (PEC) (25º16' -25º34'S; 48º17' -48º42'W), north of the state of Paraná and southern Brazil (Fig. 1), connected to the open sea by three tidal channels (Mantovanelli, 2004). The region is a considered Natural World Heritage site and a Biosphere Reserve by the United Nations since 1991 (Lana et al., 2001), and it is composed of five main bays: Paranaguá, Antonina, Laranjeiras, Pinheiros, and Guaraqueçaba. This area contains several marine protected areas, as is the case of the Ilha do Mel State Park, the Ilha do Mel Ecological Station or the Superagui National Park; the latter is considered a Natural Heritage Site, a Biosphere Reserve, and a Natural and Historical heritage of Paraná.
The PEC is one of the largest estuarine systems in the Southwest Atlantic, encompassing a total area of 612 km 2 and extensive tidal flats (295 km 2 ), with a predominant mangrove cover (Faraco et al., 2010). The average depth inside the PEC is 5.4 m and the residence time of its waters is 3.49 days (Lana et al., 2001). This estuarine system has a moderate vertical gradient of salinity and semidiurnal tides exhibiting diurnal inequality (maximum variation, 2.7 m) and consistent seasonal circulation and stratification (Marone et al., 2005). It is substantially influenced by freshwater flows from inland drainage, especially in the summer, which includes the rainy periods when the magnitude of river discharge is about five times greater than in the winter months (Mantovanelli et al., 2004).
Data collection and sample processing. Eight sampling sites oriented in the northsouth direction of the PEC were set up for this study (Fig. 1). Monthly collections took place from May 2000 to April 2001, totaling 94 samples (we could not sample in December 2000 and April 2001 at Site 2). For fish captures, we used a beach seine measuring 30 m long, 3 m high, and 5 mm stretched mesh opening. For each sampling site in each month, two replicate deployments were made at water depths between 30 cm and 1.10 m. All deployments were 50 m long and were performed parallel to the shore, with replicates separated by 5 m. Fish caught were packed in plastic bags, identified, and preserved in ice for later procedures. The collected material was identified in the laboratory using identification keys (Figueiredo, Menezes, 1978, 1980, 2000Menezes, Figueiredo, 1980, 1985. Water samples were taken using a Van Dorn bottle to measure salinity, dissolved oxygen, and temperature. Concentrations of dissolved oxygen were determined by the Winkler method, according to Grassoff et al. (1983); temperature data were obtained using a mercury thermometer and salinity was measured with a refractometer. Transparency was measured using a Secchi disk. The specimens collected for the present work were deposited in the fish collection of Museu de História Natural do Capão da Imbuia (MHNCI), Curitiba, Paraná State, Brazil (Tab. S1).
The spatial-temporal variability of the environmental predictor variables (temperature, salinity, dissolved oxygen, historical rainfall and transparency) was shown graphically. As there is no spatial variability for rainfall and spatial thermal variation was overall low, graphics with their values averaged by months were shown.
To assess the representativeness of the sampling of the fish community, a cumulative species curve was created in the vegan R package using the specaccum function (Oksanen, 2019), based on all samples collected. A modeled curve was drawn based on the Coleman's species richness estimator (Coleman et al., 1982).
Generalized Linear Model (GLM) represent an alternative to Gaussian linear models, especially when assumptions of the Gaussian models (normality of residuals, homogeneity of variance of residuals, independence) are not met (Zuur et al., 2009). GLM was used to determine the relationship between the species richness (total number of species) and the following potential predictors: time (expressed as the total number of days from the first -May 1 st , 2000 -to the last-April 30, 2001 -sampling day), distance from each sampling site to the mouth of the estuary (D, see Fig. 1 to identify the estuarine mouth), salinity, water temperature, dissolved oxygen, transparency, and 6/20 ni.bio.br | scielo.br/ni Estuarine fishes and environmental variables historical rainfall. Data referring to historical rainfall (monthly average between 1975 and 2015) were extracted from the database of the Paraná Agronomic Institute (IAPAR, 2022). Pair-plots were constructed to verify the relation among the response variables and predictors. Distance showed an approximately quadratic relation with richness, so a quadratic term of this variable was included as a potential candidate predictor (Zuur et al., 2009). A negative binomial distribution was considered using the glm.nb function of the MASS R package (Venables, Ripley, 2002).
For abundance data (total number of individuals), GLM's with both Poisson and Negative Binomial result in incorrectly specified, over disperse (see below for model validation) and low-r 2 models. Following Ives (2015)'s approach, we used a Gaussian linear model (function lm) based on log-transformed data of abundance to model the relationship among abundance and the same set of predictors used in the model richness.
For both models (richness and abundance), uniformity (KS test -to detect heteroscedasticity), overdispersion (dispersion test), outlier detection (outlier test), residuals plots were used to assess the validation of the models with the function plotSimulatedResiduals of the DHARMa R package (Hartig, 2017). With this package, temporal and spatial autocorrelation were assessed with the functions testSpatialAutocorrelation and testTemporalAutocorrelation, respectively (Hartig, 2017). For the richness model, normality of residuals was check with the Anderson-Darling test using the function ols_test_normality of the olss R package (Hebbali, 2020).
Model selection process was done using both stepwise algorithms (backward and forward) and both Akaike Information Criterion (AIC) and Bayesian Information Criterion (BIC). A total of four models emerged from this process (backward-AIC, forward-AIC, backward-BIC, forward-BIC), each one having the highest explanatory power (i.e., lowest value of the information criterion-AIC or BIC). Among the four model candidates, the best one was that with the highest value of r 2 (Venables, Ripley, 2002), calculated with the function rsq of the package R rsq (Zhang, 2022).
Prior to model selection, multicollinearity among candidate predictors was investigated run the function VIF (variation inflation factor) of the car R package (Zuur et al., 2009;Fox, Weisberg, 2019). Variables with high VIF (> 5) were excluded from the full model (Zuur et al., 2009). To show the effect of the model predictors, graphics were constructed using the function allEffects of the effects R package (Fox, 2003).
The relationship between fish assemblage structure (defined by species composition and their numerical abundances) and potential predictors (the same set used in the linear models) was assessed by using the distance-based, non-parametric linear modeling available in DistLM (Anderson et al., 2008) in the PERMANOVA+ add-on for the software PRIMER 7, v. 7.0.21. This routine fits predictor variables (matrix of predictors) to a set of response variables (the response abundance similarity matrix) on the basis of any resemblance measure. The abundance similarity matrix was built using the Bray-Curtis dissimilarity index, based on square-root transformed data of numerical abundance. The corrected Akaike information criterion (AICc) and step-wise were, respectively, the selection criterion and selection procedure used. The significance of test results was based on 9,999 permutations (Anderson et al., 2008).
The influence of the predictors of the selected model on the fish assemblage structure was visually assessed by the constrained ordination distance-based redundancy analysis (dbRDA). Vectors on dbRDA ordination determine the strength and direction of the 7/20 ni.bio.br | scielo.br/ni Luís Henrique M. C. Vergès, Riguel F. Contente, Camila Marion, Cívil P. C. del Castillo, Henry L. Spach et al. relationship between the predictors and the ordination axes. The length and direction of the vectors indicate the intensity and direction of the influence of the preditor in structuring the assemblage. To explore the relationship between species and predictors, both species and predictor variable vectors were shown on the ordination (Anderson et al., 2008). Only species showing a Person correlation coefficient (|r|) ≥0.3 with one or more axes were displayed. The R version 4.1.2 was used to run functions in R packages (R Development Core Team, 2018).

RESULTS
In general, salinity, transparency, and dissolved oxygen increased with decreasing the distance from the mouth of the estuary. While there is no clear seasonal tendency in transparency and dissolved oxygen, salinity was clearly lower in the rainy seasons than in dry ones across the estuary (Fig. 2).
Rainfall in dry seasons (April to September) was lower (< 200 mm) than in rainy seasons (November to March; > 200 mm). Water temperature was the lowest in most part of dry seasons (May-September; < 20 ºC), and highest in rainy seasons (October-March; Fig. 3). In total, 84,601 specimens of 96 taxa were caught in the shallow areas of the northsouth axis of PEC, with a predominance of Mugil sp., Atherinella brasiliensis (Quoy & Gaimard, 1825), Lycengraulis grossidens (Agassiz, 1829), Harengula clupeola (Cuvier, 1829), and Anchoa tricolor (Spix & Agassiz, 1829), which corresponded to 76.5% of the total catch. The majority of taxa (85 species) occurred with a relative frequency in the samples of less than 1% (Tab. 1). A total of 35 families were registered. The families with the highest number of species were Carangidae (13 species), Gobiidae (7), Engraulidae (6), Mugilidae (6), and Sciaenidae (6). Although there was no clear stabilization, the species accumulation curve reveals that from the 20th sample there was a gradual decrease in the slope and the among-sample variability, tending to stabilize (Fig. 4).
A January sample from site 7 at salinity 25 was disproportionately abundant (contained ~ 4,000 H. clupeola and ~ 30,000 Mugil juveniles), clearly representing an outlier; thus, it was removed from analyses. Four samples taken in salinity between 20 and 30 were also highly abundant (> 3,000 individuals) and dominated by one species: Anchoa parva (Meek & Hildebrand, 1923) (March), A. brasiliensis (August), L. grossidens (January and February). However, diagnostic analyses (see below and Fig. S2) did not identify them as outliers; thus, they were maintained for analyses. The variation inflation factor revealed that D and D 2 showed high collinearity (VIF > 5) and, therefore, both were excluded from the full richness model. The most parsimonious model displaying the highest r 2 (r 2 = 0.42) related richness to temperature, transparency, salinity, and time (Tab. 2). The model indicates that the richness tended to increase with increasing salinity, temperature and time, and tended to decrease with increasing transparency (Fig. 5).

Estuarine fishes and environmental variables
The most parsimonious model to relate assemblage structure and predictors was formed by five variables: D, temperature, transparency, salinity, and time (Tab. 3). These variables explained 26% of total variation in fish assemblage structure. The visualization of this model (given by the constrained dbRDA ordination) (Fig. 7) indicates two clear gradients of influence: D-salinity (inversely related) and temperature-time (directly related). The vector of groupings of species (Fig. 7) (Holbrook, 1847), and M. americanus Linnaeus, 1758, were more abundant in more saline lower areas; and (2) an increase in abundance of all species toward warmer waters and warmer and rainier months.

DISCUSSION
The results of the statistical models of this study support our hypothesis that richness and abundance of shallow water fish assemblage of the north-south axis of the PEC respond to (I) seasonal variability of variables related to estuarine productivity (time, temperature, and turbidity -the opposite of transparency); and (II) salinity gradient. Toward warm and rainy summer seasons, both temperature and turbidity rise, this latter because of an increase runoff in the order of 170% annual average (Lana et al., 2001). This higher land-based runoff means higher continental-source nutrient concentrations that, combined with warmer waters, boost the biological productivity, increasing the abundance and availability of prey for fishes within the estuary (Cabral et al., 2007;França et al., 2009). Many species have critical events of their life cycle (like reproduction and recruitment) tightly coupled to conditions of increased food supply in the habitat, features that maximize fish juvenile growth and survival (Robins et al., 2006). The positive relationship between such predictors/proxies of estuarine productivity and total abundance and richness should reflect the increase of use of the local by individuals and species to feed, reproduce, shelter, and settle (Elliott, Hemingway, 2002). Most of the increase in fish abundance is likely a consequence of increased rates of multi-specific recruitments (Pichler et al., 2017). The association between high productivity in the estuary and high richness and abundance has also been widely reported in other tropical (Nagelkerken, 2009) and subtropical (Harrison, Whitfield, 1995) estuaries.
According to our models, fish abundance tended to increase toward inner reaches of the estuary. This pattern may be a generalized consequence of the inner shallow areas being more favorable to occupation by small fishes (especially juveniles that need to maximize growth and survival) because, in relation to the outer reaches of the estuary, inner shallow areas have greater food supply, more shelters (mangrove roots, tidal flats) are less turbulent, and are less hydrodynamic (Vasconcelos et al., 2009).
Salinity is an important factor in structuring fish communities in estuaries (Barletta et al., 2005;Sosa-Lopez et al., 2007), defining preferred areas of species' occurrence. In our study, richness responded to salinity, increasing from the upper to lower reaches of the estuary. In general, there is more marine species using the estuaries than freshwater species (Nicolas et al., 2010;Pichler et al., 2015;Potter et al., 2015). Marine species have enhanced physiological mechanism to osmoregulation compared to freshwater species. Thus, near the ocean, high salinity waters of the lower estuary are closer to the marine pool of species, determining higher richness in these areas, thus justifying the observed pattern (Sosa-Lopez et al., 2007).
The abundance model also indicates that abundance increased with increasing salinity. Four highly abundant, one-species samples (> 3000 individuals) taking in salinity between 20 and 30 may have determined such a tendency that diverged from previous reports that found the opposite (increasing the abundance with decreasing salinity, especially in intermediary values; e.g., Bot Neto et al., 2018). Although the model coefficient of temperature, time and transparency suggested higher total abundance during the rainy season, the year-round existence of very large schools of estuarine species in higher salinity is quite common situation within PEC (e.g., Ignácio, Spach, 2009) and in other Neotropical coastal systems (e.g., Contente et al., 2020). The model coefficient for salinity may have captured this phenomenon.
The results of this study also demonstrate that the structure (composition and abundance) of the ichthyofauna of the north-south axis of the PEC was organized by the temporal, thermal, saline and spatial gradients, supporting previous studies in other estuaries Vasconcellos et al., 2010). Eucinostomus argenteus, S. testudineus, S. greeleyi, A. lineatus, and B. soporator were more abundant in less saline inner areas, whereas T. falcatus, C. faber, T. carolinus, M. littoralis, M. americanus, T. goodie, and T. marginatus were more abundant in more saline lower areas. Such patterns were also observed in other PEC environments (Hackradt et al., 2010;Contente et al., 2011) as well as in other estuaries (e.g., Babitonga Bay; Vilar et al., 2011). Matched with the results of the abundance model, these species had moderate-to-strong association with temperature and time, indicating an increase in their abundance with the proximity of the most productive seasons. The increase in the number of individuals of these species is related to the increase in their recruitment and reproductive rates since most of these species commonly occur as juvenile in the PEC marginal habitats (see Pichler et al., 2017).
It is important to note that linear models provide predictions of the response variable based on the combined effect of the predictors, and not on their isolated effects (Zuur et al., 2009). The same is true for the multivariate model (Anderson et al., 2008). This demonstrates the joint effect of the structuring factors and thus the complexity of the processes that structure the estuarine fish faunas (Hackradt et al., 2010;Pichler et al., 2017). Such a complexity was partially captured by our data, as the coefficients of determination for all models were relatively low (0.26 < r 2 < 0.42), indicating that other factors, not considered here, may play a role in structuring the fish fauna (see below).
In this study, we demonstrate the importance of the salinity gradient and variables that, directly or indirectly, play key role on the productivity cycles of the estuary to maintain the local fish biodiversity structure. The environment during the rainy period was found to be critical for the species' life cycle; therefore, local human-interventions, like dredging, should be avoided during this period (Possatto et al., 2016). Eventual shortage in freshwater and nutrient inputs due to impacts on the watershed that supplies the north-south axis of the PEC may increase salinity and reduce the local productivity (Sosa-Lopez et al., 2007). Such changes and warmer waters are also expected for the region with climate change (Pörtner et al., 2021). Salinization, low productivity, and higher water temperature are expected to affect the integrity of the local fish assemblage.
A monitoring program of the fish fauna and its structuring variables would be indispensable to verify the temporal nekton integrity. The monitoring could include variables that were not considered here but may play an important role in driving the ichthyofauna, such as prey availability, marginal structural complexity (mangroves and saltmarshes), turbulence-current, reproductive activity, slope of intertidal habitat and sediment (Martinho et al., 2009). Social actors that use (fisheries), study (local university), or negatively affect (port) the local fish biodiversity could act collectively to carry out a participatory monitoring program, which has proven to be very efficient in generating robust and useful results for conservation (Roque et al., 2018). In this line, for example, while the port authority could finance the monitoring, researchers could contribute with technical expertise and fishermen (due to their accurate ethno-ecological knowledge), with the spatial localization of the habitats that hold the greatest richness and abundance of juveniles, which should be prioritized in conservation actions.