Larval fish assemblages in nearshore waters of southeast Gulf of California: vertical and temporal patterns

The influence of salinity and temperature on larval fish assemblages, as well as, temporal and vertical patterns on larval fish assemblages off the inlet of the Presidio River, southeast coast of the Gulf of California were evaluated. Samplings for ichthyoplankton and environmental parameter measurements were carried out at three depths, in five sampling sites, during September and December 1994, and April and June 1995. Anchovies and herrings (Clupeiformes) were the most abundant larval fish accounting for 77% of the total abundance. A significant relationship between salinity and the abundance of larval herrings (Opisthonema medirastre and O. dovii) and between the water temperature and the abundance of the weakfish (Cynoscion reticulatus) was found. Anchovies (Anchoa lucida, A. walkeri, and A. nasus) were present in all sampled months, while O. medirastre and O. dovii occurred in December and June, and C. reticulatus in September. Larvae of pelagic fish were more abundant near the surface, while larvae of demersal fish were more abundant near the bottom. The present study, by emphasising the importance of considering temporal and vertical changes in larval fish assemblages in coastal environments with ecological and economic importance, will be useful for designing more efficient sampling programs.


INTRODUCTION
In coastal ecosystems, ichthyoplankton research has particular importance in ecological studies, allowing knowledge of the population dynamics of commercial species (Hale et al., 2018); determining spawning areas and times (Beck et al., 2003;Macedo-Soares et al., 2009), and taxonomy, helping to complete the life cycle of the species (Aceves-Medina et al., 2003). Understanding how larval fish communities are impacted and modulated by the physical environment they are living in is one of the key goals in marine ecology research (Chambers, Trippel, 1997), and is essential for evaluating and achieving sustainable utilisation of fishery resources.
Coastal ecosystems exhibit spatiotemporal gradients, in both environmental factors and ichthyoplankton assemblages, because of their tight physical-biological coupling (Marques et al., 2006). Environmental conditions fluctuate at different time and space scales inducing high variability on the structure and distribution of fish larvae and eggs, as well as other planktonic populations (Primo et al., 2011;Bruno et al., 2018;Zhang et al., 2019). Mixing processes such as those occurring in estuarine areas, where freshwater courses discharge into the sea, lead to a considerable vertical variation in factors like light, currents, hydrostatic pressure, and food availability (Gibson, 2003;Miller, 2007;Hale et al., 2018). These vertical changes affect the abundance, composition, and vertical distribution of larval fish (Olivar et al., 2018;Reglero et al., 2018) and consequently can affect essential ecological processes such as feeding, transport, growth, recruitment, and survival (Landaeta et al., 2015;Hale et al., 2018;Gray et al., 2019;White et al., 2019).
The Gulf of California is a region with high biodiversity along the Eastern Pacific (Brusca, 2010) where important oceanographic processes take place (Lavín, Marinone, 2003) leading to physical fluctuations and high spatial and temporal variations in abundance and species composition (Aceves-Medina et al., 2004;Urias-Leyva et al., 2018). The coast of the state of Sinaloa, Mexico, in the entrance of the Gulf of California has numerous estuarine systems of ecological and economic importance along the coastline that are subject to diverse anthropogenic pressures (overfishing, climate change, loss of biodiversity and habitats, and pollution) (Amezcua, 2014;Páez-Osuna et al., 2017). Many species of fish and invertebrates use these estuaries as spawning, feeding, and protected areas for larval and juvenile stages (Ramírez-Rodríguez et al., 2014); however, little research has been done about the composition and abundance of larval fish in these ecosystems (Álvarez-Cadena et al., 1984;Sánchez-González et al., 2018). Given the scarce information that currently exists on the ecology of ichthyoplankton in the eastern coast of the Gulf of California the goal of this work was to study the temporal and vertical variability in the larval fish community off the inlet of the Presidio river, which belongs to the Huizache-Caimanero estuarine system in southern Sinaloa. For this purpose the aims were to: i) analyse the larval fish assemblage structure considering the influence of salinity and temperature, and ii) evaluate temporal and vertical changes in the larval fish assemblage. The working hypothesis is that the larval fish assemblage includes tropical, subtropical, estuarine, and coastal species, with structural changes related to variations in salinity, temperature, sampled month, and depth.

MATERIAL AND METHODS
Sampling sites. The estuarine system of Huizache-Caimanero is classified as an intermittent estuary (Amezcua et al., 2019), which supports small-scale fisheries mainly for shrimp, but also for fin fish. The study area was located in nearshore waters in front of the Presidio River inlet (23°06'N, 106°18'W), which is part of the Huizache-Caimanero estuarine system (Fig. 1). This estuary covers an area of 175 km 2 on average, reduced to 65 km 2 during the dry season, and has two main lagoons, placed between two rivers and separated from the Pacific Ocean by a barrier island (Amezcua, 2014). Both Presidio and Baluarte rivers are connected to the lagoons by narrow and winding channels permitting the transport of freshwater into the lagoons in the wet season (June-October) and the entrance of marine water as in a typical estuary (Zetina-Rejon et al., 2003).
Ichthyoplankton samples collected during four sampling cruises on board the R/V El Puma were analysed (September 1994-autumn, December 1994-winter, April 1995spring, and June 1995. Samples were taken during the daytime from five sites along a perpendicular transect to the coastline. The first and last sites along the sampling transect were located at the following geographical coordinates (latitude -longitude): 23°4'48.71"N -106°18'21.92"W and 23°2'58.10" -106°19'56.53"W ( Fig. 1). However, sampling could not be carried out in the site closest to the coast during April and September due to technical inconveniences. At each site, ichthyoplankton samples were taken from three depths: 05 m (surface), 10 m (middle water), and 20 m (bottom).
Ichthyoplankton trawls were carried out with three Bongo nets of conical shape with an approximate length of 2.0 m, a diameter at the mouth of the net of 0.5 m, and a mesh size of 500 μm. Each net was equipped with General Oceanics digital flowmeters, previously calibrated. To prevent the bottom and half-water samples from mixing in their upward path a locking mechanism operated by a messenger was used. Each trawl lasted 10 minutes at a speed of 2.5 knots. Subsequently, the samples were placed in 4% formaldehyde for preservation until they were checked in the laboratory. In addition, in each of the sampling sites, the values of temperature and salinity were recorded through the water column by means of a General Oceanics ® Conductivity-Temperature-Depth device (CTD). A data set of monthly salinity values from October 2011 to August 2014 was also used to compare with those measured during the sampling period and to examine a current annual trend in salinity. The material was deposited in the Ichthyological Collection of the Institute of Marine Science and Limnology of the University (https://datosabiertos.unam.mx/biodiversidad).
In the laboratory, the samples were filtered and divided into 16 aliquots, using a plankton separator Folsom brand KAHLSICO, in order to estimate the abundance of each taxon (Boltovskoy, 1981). The abundance was estimated from the counts of 16 percent of each sample. Fish larvae of an aliquot of each sample were sorted and identified under a stereoscopic microscope to the lowest possible taxon, using field guides or diagnostic keys (Gasca, Suárez, 1996) and abundance was expressed as individuals per cubic meter.

Statistical analyses.
We performed descriptive statistics and profile plots on the salinity and temperature, and evaluated temporal and vertical differences in both variables using the non-parametric Kruskal-Wallis test. When the test was significant, we applied the Dunn's test for pairwise comparisons (MacFarland, Yates, 2016).
The Hill's diversity number was used to estimate the diversity (Jost, 2006), because it is more likely to have a linear relationship with environmental variables (Borcard et al., 2018). This index was used to evaluate changes in diversity and evenness related to sampled months, depths, salinity, and temperature using GLMs. The index N 1 is preferable instead of Shannon's entropy H' for interpretation through linear models and was estimated as: where H' is the Shannon entropy index: and pi is the proportion of each taxon in the sample was used to estimate taxonomic diversity.
The evenness Hill's ratio was estimated as: where N 0 is the species richness.
Temporal and vertical changes in diversity were evaluated applying generalised linear models GLMs with a Gamma error distribution structure and a log-link function that is appropriate for continuous response variables and when the residuals of the model are expected to have a positively skewed distribution (Fox et al., 2015). Assumptions were checked by plotting the residuals versus the fits. Univariate GLMs were applied in the R software using the function "glm" from the Stats package in R (R Core Team, 2019), and applying a Likelihood Ratio test (LR).
Statistical analyses were also conducted to test shifts in the larval fish assemblage related to salinity, temperature, depth, and sampled month taking into account a distance-based approach using permutational multivariate ANOVA (PERMANOVA) (Anderson, 2001) and a canonical db-RDA analysis to graphically show the results (Borcard et al., 2018). In addition, a model-based approach (Warton et al., 2015) using multivariate and univariate GLMs  was applied. This analysis allows for community-level (multivariate) and species-level (univariate) inference using resampling based hypothesis testing in a more robust statistical framework .
The db-RDA was performed using Bray-Curtis distances on fourth-root transformed data to enhance the contribution of the less abundant taxa (Clarke, Warwick, 2001). The water temperature (°C), and salinity were the explanatory variables, whereas the months (September, December, April, and June), and depth (5 m, 10 m, 20 m) were the fixed factors. The test was applied using the function "capscale" from the R package "vegan", and the global significance was tested by permutation with the function "anova.cca" from the same package (Oksanen, 2019). In addition, we use a two-way PERMANOVA (Anderson, 2001) to examine the effects of the temporal and vertical variations (fixed factors) on the ichthyoplankton community. Terms found to be significant were examined individually using pairwise tests of the PERMANOVA test. As this analysis can sometimes be affected by differences in within-group variation, the data were also tested for homogeneity of group dispersions (Anderson et al., 2008). We performed the PERMANOVA in the Primer-e® software using a "Type III" sum of squares, which is appropriate for unbalanced models (Legendre, Legendre, 2012).
As a complement to the distance-based approach, multivariate and univariate GLMs were applied to test the effect of the temporal and vertical variations, salinity, and temperature on larval fish assemblage, using the statistical package mvabund in R . GLMs were fitted under the assumption of a negative binomial distribution of errors and applying a logit link function (Hilbe, 2011). Assumptions were checked by plotting the residuals versus the fits .

Temporal and vertical variations in salinity and temperature.
In the study area, water temperature (mean ± SD = 25.8 ± 4.0 °C) ranged from 19.2 to 31.0 °C. The water column was significantly warmer in September and June than in December and April (Dunn's test p-value < 0.05, Fig. 2a). Salinity (mean ± SD = 32.9 ± 3.6) ranged from 23.8 to 36.7 and was significantly highest in April and lowest in September (Fig. 2a). In general, the water temperature was homogeneous in the water column, except in June when some changes were observed both among depths and sampling sites (mean ± SD = 26.5 ± 2.4, range = 23.1-29.0); while salinity showed more spatial fluctuations in September (mean ± SD = 26.9 ± 2.4, range = 23.8-32.3) (Fig 3). However, there were no significant differences in salinity and temperature among depths (p-values > 0.05). The current annual trend in salinity showed a pattern with the highest values in May and the lowest in September (Fig 2b). The salinity average in the whole current period was 33.0 ± 2.7 (mean ± SD).
Taxonomic composition and diversity. The larval fish assemblage included tropical, subtropical, estuarine, and coastal species. A total of 17 fish species belonging to 10 families (Engraulidae, Clupeidae, Gerridae, Sciaenidae, Carangidae, Albulidae, Serranidae, Ephipidae, Pleuronectidae, and Paralichthyidae) and 4 orders (Clupeiformes, Perciformes, Albuliformes, and Pleuronectiformes) were collected. The order Clupeiformes accounted for 77% of total abundance and was represented by five species, including herrings and anchovies; whereas the order Perciformes accounted for only 21% of total abundance and was represented by ten species. Anchoa lucida (Jordan, Gilbert, 1882) was the most abundant species with an average abundance value of 3.0 ± 0.8 ind. m -3 and a maximum of 40.0 ind. m -3 . The mean, maximum, and standard error values for all species identified are summarised in Tab. 1.
The diversity Shannon index (H') had a mean value of 1.3 ± 0.5 (mean ± SD) and ranged from 0.0 to 2.2. Diversity (Hill's diversity number N 1 ) was found to be significantly greater in June than in April (p-value = 0.04), while evenness (Hill's evenness ratio E 1 ) was higher in April than in the remaining months (Fig. 4). Diversity and evenness did not vary significantly along salinity, temperature, and depth gradients.
Effects of temperature and salinity on larval fish assemblages. There was a strong FIGURE 2 | A. Box plots showing temporal variations in temperature and salinity during the sampling period (1994)(1995)  influence of salinity and temperature variations on the ichthyoplankton composition and abundance (Fig. 5). This figure is a db-RDA correlation triplot which reflects the correlations between species and continuous explanatory variables (salinity and temperature); this relationship was highly significant (p-value = 0.001) after 999 permutations. In the db-RDA analysis 48% of the total variance (the fitted variance) was explained by the explanatory variables, and 78% of this proportion is represented by the first two canonical axes, which are linear combinations of the explanatory variables. Salinity was correlated to the abundance of Opisthonema medirastre Berry, Barrett, 1963 andO. dovii (Günther, 1868), which were placed close to samples taken during June in the ordination. Temperature was correlated to the abundance of Cynoscion reticulatus (Günther, 1864), Menticirrhus elongatus (Günther, 1864), Larimus effulgens Gilbert, 1898, and Umbrina analis Günther, 1868 (Fig.  5). Based on the community-level analysis using the multivariate GLM significant effects of salinity (LR: 28.89, p-value: 0.04), but not temperature (LR: 21.99, p-value: 0.16) were observed on the whole larval fish community. Regarding the species-level analysis using the univariate GLMs the abundance of both Opisthonema medirastre and O. dovii increased with higher salinity values (LR = 8.94, p = 0.05 and LR = 9.65, p = 0.03 respectively), while the abundance of the remaining species did not fluctuate significantly with changes in temperature and salinity. Temporal patterns in the larval fish composition. Based on the community-level analysis, both the two-way PERMANOVA (Tab. 2) and the multivariate GLM (LR = 186.5, p-value = 0.001) revealed significant effects of the temporal variation on the structure of larval fish assemblage. No significant effects of the interaction between the sampled month and depth on the larval fish assemblage were found (Tab. 2).
A major temporal pattern in the abundance of both O. medirastre and O. dovii was observed (Fig. 6). A significant effect (p < 0.001) of the temporal variation on the abundance of both species was found by the univariate GLMs; these species were present only in June and December, and their mean abundance peaked in June (Fig. 6). A temporal pattern was also observed in the abundance of A. lucida; A. walkeri Baldwin, Chang, 1970;A. nasus (Kner, Steindachner, 1867); and Cynoscion reticulatus. The mean abundance of larval anchovies reached its highest values (11.00, 6.60, and 4.40 ind.m -3 , respectively) during September near the surface, while the mean abundance of larval C. reticulatus peaked in June near the bottom (2.14 ind.m -3 , Fig.6). However, these temporal patterns were not significant based on the species-level analysis using the univariate GLMs (p > 0.05).

Vertical distribution of larval fish.
In the community-level analysis PERMANOVA (Tab. 2) revealed differences in the larval fish assemblage among depths. Results yielded by the pair-wise test of the PERMANOVA showed that the community at 20 m was different from that at 5 m (t = 2.01, p-value = 0.018) and 10 m (t = 1.7, p-value = 0.049). However, there were no significant differences in the community structure between 5 m and 10 m. The variations in the larval fish assemblage across depths can be seen in Fig. 7. All three species of anchovies (A. lucida, A. walkeri, and A. nasus) were more abundant in the upper parts of the water column while larval sciaenids (C. reticulatus, M. elongatus, U. analis, and L. effulgens) were more abundant in deeper parts. However, FIGURE 7 | Vertical distribution of the abundance of the species found in the study. Bars indicate the mean value (individuals m -3 ) and error bars the Standard Error. no significant vertical variation in abundance was observed for all species, based on the univariate GLMs.

DISCUSSION
Although there are some studies on ichthyoplankton assemblages in the Gulf of California (Aceves-Medina et al., 2003, 2004, none of them have been undertaken on its southern coast. Much of the research on coastal ichthyology in the southeastern coast of the Gulf of California has focused mainly on adult fish populations (Amezcua, 2014). Our current findings expand and complement the knowledge in the ichthyology of the zone by an intensive study of the structure of larval fish communities and their vertical and temporal dynamic. Temporal variability in the diversity of larval fish could be related to the circulation in the entrance area of the gulf (Färber-Lorda et al., 2010;Urias-Leyva et al., 2018). Circulation studies at the mouth of the gulf indicate a complex flow pattern, with spatial and temporal variations in the inflow and outflow (Castro et al., 2006;Lavín et al., 2014). It has been suggested that the inflow occurs next to Sinaloa at the east side and that outflow takes place along Baja California at the west side (Lavín, Marinone, 2003). In the present work diversity generally reached the highest values in the period of the year when tropical waters flow northward to the entrance of the gulf. This flow is caused by the effect of a branch of the North Equatorial Countercurrent (Lavín et al., 2014) and transports tropical water masses with a greater diversity of planktonic species (Gasca, Suárez, 1996;Franco-Gordo et al., 2003) leading to an increased diversity in the studied area. In contrast, the lowest values of diversity match with cold-water masses moving southward from higher latitudes (Franco-Gordo et al., 2003;Lavín, Marinone, 2003). These findings are in general agreement with those obtained by Franco-Gordo et al. (2003) in the Pacific coast of Jalisco (southward Sinaloa) that determined a temporal diversity pattern over the period 1995-1996 with three main periods, one from January to May representing the influence of the California Current, a second period (July to November) representing the tropical affinity season, and a third one (June and December months) representing the transition months.
In the present study the main taxa in terms of abundance were pelagic species such as herrings (Clupeidae) and anchovies (Engraulidae), while demersal fishes in general were found to be less abundant. Species belonging to these groups are common along the Sinaloa coast (Ramírez-Rodríguez et al., 2014) and other tropical and subtropical coastal environments (Tzeng et al., 1997;Brusca, 2010;Vega-Corrales, 2010). All three species of anchovies (A. lucida, A. walkeri, and A. nasus) occurred in all sampled months while herrings (O. medirastre and O. dovii) were found in December (temperate and dry season) and June (warm and dry season) when they reached the highest abundances. Temporal changes in abundance may reflect differences in the abundance and spawning activities of adult populations and larval survival (Espinosa-Fuentes, Flores-Coto, 2004;Ramos et al., 2006;Sabatés et al., 2007). In the Pacific coast of Costa Rica, it has been reported that the discontinuity in the catch of the O. medirastre throughout the year could be related to temporal variations in water temperature and phytoplankton abundance (Funes-Rodríguez et al., 2001;Vega-Corrales, 2010).
The relationship between the distribution of planktonic organisms and environmental variables has been studied in coastal environments around the world. Salinity and temperature are frequently the most important parameters affecting the distribution and abundance of organisms (Paul et al., 2016;Santos et al., 2017;Reglero et al., 2018). In the present work salinity variations were determined by the temporal variation (Fig.  2), while no significant variations between depths were observed. The influence of the Presidio river discharge on the spatial salinity distribution was evident only during the wet warm-season (September, see Fig. 3); however, a narrow gradient of salinity (from 25 to 36) was observed along the transect. Salinity is a key physical factor that constrains the spatial distribution of estuarine ichthyoplankton (Santos et al., 2017); salinity variations have been found to be a major driver influencing fish assemblage structure in systems experiencing wider ranges of salinity such as estuaries, where ichthyoplankton assemblages may be composed of a mix of marine, estuarine, and freshwater taxa (Sánchez-González et al., 2018).
The interannual variability in salinity was influenced by changes in the flow of the Presidio River determined by the wet and dry seasons, and it was consistent with the trend observed during the period 2011-2014 (Fig. 2b). As both the mean salinity value and the interannual variability are consistent over time our findings may be considered important for the time of the study but also now.
Marine fish larvae, like other zooplankton, rarely distribute randomly in the water column; instead, they display patterns of vertical stratification depending on the species and environment (Hale et al., 2018). The vertical distribution of larval fish found in our study seems to be determined by the vertical distribution of adult populations (Gray, Miskiewicz, 2000;Franco-Gordo et al., 2002;Brusca, 2010). Larvae of pelagic fish (anchovies) were more abundant at the upper parts of the water column, while larvae of demersal fish were more abundant at deeper parts (Fig. 7). We did not find depth-related variations in the diversity index (p > 0.05). However, a pattern with a higher number of species occupying deeper parts of the water column has been reported for other coastal areas (Cha et al., 1994;Espinosa-Fuentes, Flores-Coto, 2004), which could be related to changes in factors such as light, hydrostatic pressure, and food availability (Fiksen et al., 2007;Miller, 2007;Hurst et al., 2009). The present study, by emphasising the importance of considering temporal and vertical changes in the larval fish assemblage in coastal environments with ecological and economic importance, will be useful for designing more efficient sampling programs.
Larval fish assemblages found in the present study were constituted for both tropical and subtropical species. The highest abundances corresponded to coastal herring species such as O. medirastre and O. dovii and anchovy species that inhabit both coastal and estuarine areas (A. lucida, A. walkeri, and A. nasus). Anchovies were more abundant in September, while herrings were more abundant in June and December. Moreover, larval fish of demersal species (C. reticulatus, M. elongatus, U. analis, and L. effulgens) were more abundant in deeper parts of the water column (20 metres depth) while anchovies were more abundant at the upper parts (5 and 10 metres depth). According to our results, the water temperature was positively correlated to the abundance of sciaenidae demersal species (C. reticulatus, M. elongatus, U. analis, and L. effulgens) while salinity was positively correlated to the abundance of herrings (O. medirastre and O. dovii). Therefore, our results suggest that the abundance of these species would be affected by changes in salinity and temperature in a global warming context. These kinds of studies can help better understand potential changes in the distribution of fish species considering potential changes in salinity and current systems due to events such as climate change.

ACKNOWLEDGMENTS
We thank Instituto de Ciencias del Mar y Limnología of the Universidad Nacional Autónoma de México for the funding to undertake this project.