INTRODUCTION
Lagoon-estuarine systems are formed by a complex matrix of intertidal habitats with distinctive environ-mental features that play an essential role in support of high fish production (McLusky & Elliott, 2004; Yáñez-Arancibia et al., 2004; França et al., 2012). In these systems, geomorphological processes originate geoforms (e.g., lagoons, mouths and channels), which are habitats with different physical attributes (e.g., depth, width, substrate) that influence the patterns of composition, richness, abundance, migration and the size-specific spatial use by organisms (Bartholomew et al., 2000; Saintilan, 2004; Franco et al., 2006; França et al., 2009; Hicks et al., 2010). For example, the mouth and its tidal channels are ecotones connecting the marine with the brackish environment (Webster, 2011; Lacerda et al., 2014). Furthermore, other geoforms such as lagoons, which are partially enclosed areas with low water circulation and more structural complexity (e.g., muddy bottoms, riverine and submerged vegetation), are reported as critical habitat for juvenile fish feeding and nurseries (Neves et al., 2013; Verdiell-Cubedo et al., 2013).
Habitat features have been integrated in estuarine fish research as a key factor to understand the structure of ecological assemblages and their relationship with highly variable spatial conditions and allows to define connectivity mechanisms between habitats and local migration patterns (Sindilariu et al., 2006; Franco et al., 2009; Neves et al., 2011), as has been shown by several studies carried out in estuaries and coastal lagoons of Europe, South Africa, Australia and South America (Patterson & Whitfield, 2000; Griffiths, 2001; Franco et al., 2006; Neves et al., 2011; Loureiro et al., 2016).
Some studies have been conducted on fish assem-blages of Mexican coasts, concerning mainly lagoon-estuarine systems from the Gulf of Mexico and the northwest and central Pacific areas (Warburton, 1978; Mendoza et al., 2009; Castillo-Rivera et al., 2010; Rodríguez-Romero et al., 2011; Ayala-Pérez et al., 2014). These studies indicate that physical and chemical changes (e.g., variations in salinity, water temperature, turbidity, dissolved oxygen, etc.) associated with seasonality constitute a significant driver of fish assemblage patterns. However, in the south Pacific of Mexico there are no studies that include the effect of habitat complexity on fish assemblages.
The Gulf of Tehuantepec in the Mexican Southern Pacific is one of the most important fishing areas (Bakun et al., 1999; Ortega-García et al., 2000; Wilkinson et al., 2009). This region has an extensive belt of lagoon-estuarine systems surrounded by mangroves that shelters one of the richest areas for fish fauna in the Tropical Mexican Pacific (Gómez-González et al., 2012; González-Acosta et al., 2017). The environmental differences between the geomorphic habitats that integrate these systems make them good models to test the spatial effect of its characteristics on biotic assemblages.
The rapid degradation and loss of coastal ecosystems in recent decades has increased interest in fish research and their habitat in estuarine environments (França et al., 2012; Sundblad & Bergström, 2014). The knowledge of geomorphic habitats and their effect on fish assemblages constitutes an important baseline for the development of monitoring programs and integrate ecosystem-based management of coastal resources (Pérez-Ruzafa et al., 2007; Franco et al., 2009; Sheaves et al., 2012). Therefore, the specific aims of this study were to i) determine the spatio-temporal variations of fish assemblages in a shallow lagoon-estuarine system and compare its structure in three geomorphic habitats (lagoon, channel, inlet); and to ii) determine their relationships with hydrological parameters (e.g., salinity, dissolved oxygen, depth) and physical attributes (substrate, vegetation, shelter). The tested hypothesis is that fish assemblages are likely to differ among these habitats, with the most complex areas supporting comparatively higher fish abundance, richness and a higher number of small-sized organisms compared to less complex areas, along spatial and seasonal gradients.
MATERIALS AND METHODS
Study area
La Joya-Buenavista lagoon-estuarine system (LJB) is located at the northeastern corner of Gulf of Tehuantepec, southwest Chiapas, Mexico (15°48’-15°56’N, 93°32’-93°47’W); it is a shallow water body with an area of 47.5 km2 (Fig.1). The LJB presents a sandy barrier and tides predominantly semidiurnal, with a range of 1.2 m (Lankford, 1977; Contreras, 2010). The water circulation into the system is reduced and depends mainly on the tides and the seasonal freshwater input provided by a few streams. The climate in this area is characterized by a dry season (November-April) and a rainy season (May-October), with a total annual rainfall of 1,441 mm and an average temperature of 28°C. The salinity can vary drastically in some areas, from oligohaline (salinity <5) at the peak of the rains to hyperhaline (salinity >40) in highly protected sites during the dry season (Contreras & Zabalegui, 1991).

Figure 1 Map of the study area in the La Joya-Buenavista lagoon-estuarine system, Mexico. Symbols represent sampling sites coded by habitat type. Circles: lagoon, squares: channel, triangles: inlet.
This study analyzed the fish assemblages of three geomorphic habitats in the LJB system: lagoon, channel and inlet (Fig. 1). These habitats were identified mainly based on the qualitative description of physical and morphological characteristics. The depth of the lagoon varies from 0.2 to 1.5 m, with a predominantly muddy substrate, with sparse submerged tree trunks and margins comprised of grass and fringe/basin mangroves (Avicennia germinans L. and Conocarpus erectus L.). The channel has a total length of 23 km and varies in width from about 570 m at the mouth to 200 m at the head. Most of the channel zone is surrounded by riverine mangroves (Rhizophora mangle L.). Substrate composition changes from muddy, with some sandy patches in the channel, to predominantly sandy in the mouth. The inlet flows into the Gulf of Tehuantepec through an estuary mouth of 157 m of width, the zone most exposed to tidal currents.
Data collection
Sampling was performed bimonthly from 9:00 to 16:00 h, at 13 sites during June 2013-April 2014. Sampling was conducted using monofilament cast nets (4 m diameter with 1.27 cm mesh). This gear could be used in all sites under its temporal variations in depth and habitat conditions, although an a priori characterization focused on small and juvenile components was assumed (Sheaves et al., 2007; Stein III et al., 2014). The cast net was operated from the shoreline or a boat by the same individual. Sampling effort was standardized to 10 cast nets. This number of deployments was determined by estimating the values of the asymptotes of plot-based species accumulation curves for each habitat in three pilot sites (R2 = 0.99, 0.94, 0.96). Cast net deployments were done in a radius of up to 50 m, thus avoiding the disturbance effect of three castings at the same point.
The fish specimens were weighed, measured, fixed in 10% formalin, washed in tap water after 48 h, and transferred to 70% ethanol. All fish were identified to species and counted. Scientific names and authorities were corroborated following Fricke et al. (2019). Voucher specimens were deposited at the ichthyo-logical collection of the Museum of Zoology of the University of Science and Arts of Chiapas (MZ-P-UNICACH) in Tuxtla Gutiérrez, Mexico. For descriptive purposes, selected dominant species were allocated to four ecological categories according to their salinity tolerance (Myers, 1966; Castro-Aguirre et al., 1999; Day Jr. et al., 2012). These species were also allocated to functional guilds according to their trophic group following the criteria of Elliott et al. (2007), with information from Froese & Pauly (2016). Four ecological groups were identified in this study: estuarine (fish that are residents in brackish waters), marine euryhaline (can tolerate a wide range of salinity), marine stenohaline (can tolerate a narrow range of salinity) and secondary freshwater (salt-tolerant inland fish). Likewise, five trophic groups were identified: piscivores, which are species whose diet consists mainly of fish; zooplanktivores, species that feed predominantly on zooplankton; zoobenthivores, species that feed on bottom invertebrates; detritivores, species that feed on the bottom, selecting fine particles, benthic diatoms, meiofauna and sediments; and omnivores, which feed on algae and a variety of invertebrates.
Hydrological variables were measured on each sampling event before the fish collection. Salinity, water temperature (°C), dissolved oxygen (mg L-1), and pH were determined at the mid-depth of each site using a multiprobe meter (YSI 556). The depth (to the nearest 0.1 m) and water clarity (Secchi's depth in cm) were also determined.
Habitat structure descriptors (substrate, vegetation and shelter) were examined in the system (in this study, habitat structure was considered as the physical attributes that may provide protection and food for fish species; Green et al., 2012). At each site, a plot of 1 m2 was randomly positioned to visually evaluate the substrate composition as the percentage covered by muddy, sandy or mixed deposits. The presence of other physical elements as tree trunks and riprap was also recorded. Values for vegetation were assigned based on the dominant marginal coverage (>60% of mangroves or grasses) in a 200 m2 area adjacent to the sampling site. The shelter was assessed according to the distance of the sites with the inlet or the channel. Following Neves et al. (2013), each habitat attribute was assigned a value from 1 to 3 to qualify structure, and these were added to achieve a total score, estimated to obtain the overall degree of physical habitat structure. Each site was categorized with ranges from 3 (lowest habitat structure) to 9 (highest habitat structure). For example, highly exposed sites (such as the mouth) with bare sand substrate and lacking plant (e.g., grasses) or physical structures were qualified with a value of 3, whereas other sites sheltered from the currents (such as the lagoon), with a combination of muddy substrate, submerged trunks and coverage dominated by mangroves would get a score of 9.
Data analysis
The species dominance was determined using the Importance Value Index (IV), estimated as a percentage of the numerical abundance, biomass and frequency of all species (Krebs, 1999). We considered species that, in total, had 75% of the IV as dominant (Peralta-Meixueiro & Vega-Cendejas, 2011).
The abundance data were fourth-root transformed to minimize the influence of overly abundant taxa. The data set was divided by the three geomorphic habitats (lagoon, nine sites, inlet and channel, two sites each; Fig. 1) and months were summarized in seasons as pre-wet (June), wet (August-October), pre-dry (December) and dry (February-April). An individual-based rarefaction method was used to estimate the expected number of species in a sample [E (Sn)] and calculate an abundance-corrected dataset to compensate for differences in sampling effort between habitats and seasons. The E (S) was generated for a constant number in a random set of m individuals from the reference sample (m < n) (Colwell et al., 2012). The variation of the rarefied species richness, as well as that of the environmental parameters between habitat and season, were analyzed using one-way analysis of variance (ANOVA), having verified the assumptions of normality and homogeneity of variances. Non-parametric ANOVA (Kruskal-Wallis test) was used in cases where heteroscedasticity was found even after the logarithmic transformation of data. A Bonferroni test (or Mann-Whitney pairwise to non-parametric test) followed the ANOVA procedures every time that the null hypothesis was rejected at α = 0.05 (Zar, 2010).
Non-parametric multivariate analyses were conducted to assess the spatio-temporal changes in fish assemblage structure (species composition and abundance) and their relationship to environmental parameters, after evaluating the homogeneity of the multivariate dispersion (PERMDISP) (Anderson, 2006). Preliminarily, a permutational multivariate analysis of variance (PERMANOVA) (Anderson, 2001), conducted on a Euclidean similarity matrix of all environmental data, assessed differences, and the interaction between habitats and seasons. Furthermore, principal component analysis (PCA) was performed on the normalized hydrological data (log x+10) to explore gradients in the sampling events. To avoid bias in the PCA related to highly correlated hydrological data, before PCA, a Pearson correlation test was performed, and variables that show high collinearity (r > 8) were not considered in this analysis. Later, differences in fish assemblages between habitat and seasons were also assessed by PERMANOVA, including habitat type (fixed factor), seasons (fixed factor), sampling sites (nested in habitat type) and months (nested in seasons), conducted on a Bray-Curtis similarity matrix. A pairwise post-hoc comparison was performed, while results off that PERMANOVA were significant.
The similarity percentage analysis (SIMPER) (Clarke, 1993) was used to describe the variance explained per species by habitat. Standard length (SL) data of dominant species selected by SIMPER were compared with a Kruskal-Wallis test and the Mann-Whitney pairwise test to assess possible variations in the size of individuals per habitat.
Finally, distance-based linear modelling (DistLM) (Anderson et al., 2008) with a step-wise procedure and an adjusted R2 selection criterion (Legendre & Anderson, 1999), was performed on the fourth-root abundance matrix and the normalized environmental data (including habitat structure values). The resulting model was represented in a multi-dimensional space with a distance-based redundancy analysis (dbRDA) (McArdle & Anderson, 2001), to assess the relative contributions of environmental variables in structuring fish assemblage and detect a possible spatial gradient. We worked with a 5% significance level (P ≤ 0.05) in all statistical tests. All multivariate analyses were performed with routines in PERMANOVA+ for PRIMER 6 software (Clarke & Gorley, 2006; Anderson et al., 2008). All other analyses were carried out with package Vegan in language R, version 3.1.2 (Oksanen et al., 2014; R Development Core Team, 2014).
RESULTS
Environmental variation
The PERMANOVA results showed that hydrological variables were significantly different between habitats (F = 9.25, P < 0.01), seasons (F = 18.98, P < 0.01) and the interaction (habitat × season, F = 4.69, P < 0.01). Likewise, differences in some hydrological parameters were consistent (P < 0.05) between seasons and habitats, or both. In paired comparisons within temporal factor, major differences between habitats were observed in all seasons (P < 0.02), except for the inlet-channel pair; oxygen saturation was low during the dry season in the channel (1 mg L-1; April), while higher concentrations were for the lagoon (5.3-9.5 mg L-1) in the wet and pre-dry seasons. According to an adaptation of the Venice system used to describe salinity zones in estuaries (Whitfield, 1998), the inlet and channel areas are typically polyhaline, ranging from 22.1 to 29.5, whereas the lagoon was significantly different from mesohaline conditions (13). Seasonal variations in salinity were recorded, ranging from oligohaline (wet season; 3.3) to marine conditions (dry; 34.1). Water temperature and pH presented seasonal variations, while the depth and water clarity only showed changes between habitats. The habitat structure score presented its highest mean value for the lagoon (8.6), while the lowest value was obtained for the inlet (3.3) (Table 1).
Table 1 Means (and standard deviation) of hydrological parameters in the La Joya-Buenavista lagoon-estuarine system (Mexico), June 2013-April 2014, and results of the ANOVA (F) and Kruskal-Wallis test (H) analyses for spatial and temporal comparisons for each parameter.
Hydrological parameters | Season | Inlet | Channel | Lagoon | Spatial | Temporal |
Water temperature (°C) | ||||||
Pre-wetᵃ | 31.4 (0.1) | 32.0 (0.4) | 33.1 (1.6) | H = 2.99 | F = 36.84 | |
Wetᵃ | 30.5 (0.1) | 30.1 (1.3) | 31.5 (2.5) | P = 0.224 | P = 0.000 | |
Pre-dryb | 28.6 (0.0) | 27.5 (0.1) | 25.8 (1.7) | |||
Dryᵃ | 30.7 (1.1) | 32.6 (1.3) | 32.3 (1.4) | |||
Depth (cm) | ||||||
Pre-wet | 89.4 (0.4) | 50.3 (4.2) | 55.3 (16.1) | H = 7.3 | H = 3.21 | |
Wet | 66.3 (8.5) | 82.7 (64.9) | 80.5 (52.5) | P = 0.026 | P = 0.361 | |
Pre-dry | 147.5 (0.0) | 181.4 (49.5) | 73.2 (28.1) | I, C > L | ||
Dry | 80.2 (77.5) | 96.1 (35.7) | 44.1(24.4) | |||
Water clarity (cm) | ||||||
Pre-wet | 87.4 (0.2) | 47.5 (4.2) | 52.4 (16.1) | H = 7.31 | H = 3.21 | |
Wet | 66.3 (8.5) | 40.0 (11.9) | 55.4 (20.0) | P = 0.026 | P = 0.361 | |
Pre-dry | 145.0 (1.8) | 125.6 (49.5) | 40.0 (18.1) | I, C > L | ||
Dry | 69.5 (77.1) | 86.3 (35.7) | 44.5 (24.0) | |||
Dissolved oxygen (mg L-1) | ||||||
Pre-wetab | 4.6 (0.1) | 2.4 (1.0) | 4.8 (1.2) | F = 25.08 | F = 4.45 | |
Wetb | 4.4 (0.0) | 4.4 (0.5) | 5.4 (1.4) | P = 0.000 | P = 0.006 | |
Pre-dryab | 3.4 (0.3) | 2.5 (0.3) | 6.5 (1.9) | L > I, C | ||
Dryac | 1.9 (0.9) | 2.4 (1.8) | 4.6 (1.2) | |||
Salinity | ||||||
Pre-wetᵃ | 31.9 (0.1) | 20.3 (1.5) | 14.1 (5.8) | F = 8.37 | F = 43.32 | |
Wetb | 23.6 (8.9) | 16.7 (5.9) | 3.3 (2.4) | P = 0.006 | P = 0.000 | |
Pre-dryab | 29.1 (0.6) | 17.5 (2.1) | 7.1 (2.3) | I, C > L | ||
Dryc | 34.0 (0.5) | 34.1 (0.4) | 27.7 (5.6) | |||
pH | ||||||
Pre-wetᵃ | 8.4 (0.1) | 8.2 (1.0) | 8.9 (1.0) | F = 2.78 | F = 10.37 | |
Wetb | 7.5 (0.5) | 7.6 (0.3) | 7.1 (0.7) | P = 0.068 | P = 0.000 | |
Pre-drybc | 7.4 (0.1) | 7.6 (1.0) | 7.7 (0.7) | |||
Dryac | 7.1 (1.2) | 6.7 (0.3) | 8.8 (0.8) | |||
Inlet | Channel | Lagoon | ||||
Substrate | 1.1 (0.1) | 2.2 (0.2) | 2.9 (0.2) | |||
Vegetation | 1.2 (0.3) | 3.0 (0.0) | 2.8 (0.5) | |||
Shelter | 1.1 (0.1) | 2.1 (0.1) | 2.9 (0.1) | |||
Mean value | 3.3 (0.4) | 7.3 (0.4) | 8.6 (0.5) |
Significant P-values (<0.05) are in bold. Seasons having the same superscript letter were not significantly different within each temporality. Habitat structure values also are presented. I: inlet; C: channel; L: lagoon.
The PCA of the hydrological variables explained 63.6% of the variance for the first two axes (Fig. 2). This model considered only four explanatory variables after reducing multicollinearity. The PC1 contributed with 36.8% of the explained variance and it was strongly related (based on eigenvectors of correlation matrix) with depth (-0.56) and salinity (-0.45); whereas, the PC2 contributed with 26.8% of the explained variance and better correlated with dissolved oxygen (0.43) and pH (-0.41).

Figure 2 Principal component analysis environmental biplot, La Joya-Buenavista lagoon-estuarine system, Mexico. The length and direction of the vectors represent the strength of the relationship concerning a unit circle. Symbols represent an individual sampling event coded by season. Gray circles: pre-wet, gray inverted triangles: wet, white circles: pre-dry, black triangles: dry.
Fish assemblage structure
Two thousand sixty-four fish individuals were collected by cast net (from 76 sampling events during the study period) and were identified in 48 species and 21 families. The fish assemblage in the LJB was dominated by Lile gracilis (13.91%), followed by Astatheros macracanthus (10.08%) and Eucinostomus currani (9.98%); only three species, Mugil curema, Dormitator latifrons and Gerres simillimus contributed with an amount as large as 36.78% of the total biomass. According to the IV, 15 species together contribute 75% of relative value (Table 2).
Table 2 Relative importance value of the fish species for the La Joya-Buenavista lagoon-estuarine system (Mexico).
Family | Specie | n | SL (cm) | %n | %B | %F | %IV |
---|---|---|---|---|---|---|---|
Clupeidae | Lile gracilis Castro-Aguirre & Vivero, 1990 | 287 | 3.3-8.2 | 13.91 | 5.52 | 6.42 | 8.61 |
Cichlidae | Astatheros macracanthus (Günther, 1864) | 208 | 2.1-13.4 | 10.08 | 9.44 | 5.88 | 8.47 |
Mugilidae | Mugil curema Valenciennes, 1863 | 129 | 4.3-21.7 | 6.25 | 12.41 | 4.28 | 7.65 |
Gerreidae | Eucinostomus currani Zahuranec, 1980 | 206 | 1.5-9.5 | 9.98 | 8.13 | 4.81 | 7.64 |
Eleotridae | Dormitator latifrons (Richardson, 1844) | 115 | 1.7-15.7 | 5.57 | 14.34 | 2.67 | 7.53 |
Gerreidae | Gerres simillimus Regan, 1907 | 82 | 1.9-18.3 | 3.97 | 10.03 | 5.35 | 6.45 |
Gerreidae | Diapterus brevirostris (Sauvage, 1879) | 113 | 2.3-8.8 | 5.47 | 5.90 | 6.42 | 5.93 |
Poeciliidae | Poecilia nelsoni (Meek, 1904) | 200 | 1.7-5.1 | 9.69 | 2.30 | 3.21 | 5.07 |
Centropomidae | Centropomus robalito Jordan & Gilbert, 1882 | 121 | 3.1-16.7 | 5.86 | 3.48 | 4.28 | 4.54 |
Cichlidae | Amphilophus trimaculatus (Günther, 1867) | 81 | 2.3-12.7 | 3.92 | 1.54 | 3.21 | 2.89 |
Atherinopsidae | Atherinella guatemalensis (Günther, 1864) | 63 | 1.6-6.5 | 3.05 | 0.15 | 4.28 | 2.49 |
Carangidae | Oligoplites altus (Günther, 1868) | 16 | 6.5-18.7 | 0.78 | 2.95 | 2.67 | 2.13 |
Poeciliidae | Poecilia sphenops Valenciennes, 1846 | 75 | 2.0-5.1 | 3.63 | 1.14 | 1.07 | 1.95 |
Engraulidae | Anchoa mundeola (Gilbert & Pierson, 1898) | 26 | 4.0-8.4 | 1.26 | 0.60 | 3.74 | 1.87 |
Gobiidae | Gobionellus microdon (Gilbert, 1892) | 23 | 5.3-15.0 | 1.11 | 1.19 | 3.21 | 1.84 |
Eleotridae | Gobiomorus maculatus (Günther, 1859) | 29 | 7.5-14.1 | 1.41 | 2.58 | 1.07 | 1.68 |
Clupeidae | Opisthonema libertate (Günther, 1867) | 47 | 5.2-7.9 | 2.28 | 1.18 | 1.07 | 1.51 |
Gerreidae | Eugerres axillaris (Günther,1864) | 9 | 7.3-14.1 | 0.44 | 1.53 | 2.14 | 1.37 |
Centropomidae | Centropomus armatus Gill, 1863 | 20 | 6.3-14.4 | 0.97 | 1.49 | 1.60 | 1.36 |
Atherinopsidae | Membras gilberti (Jordan & Bollman, 1890) | 27 | 8.3-11.1 | 1.31 | 1.49 | 1.07 | 1.29 |
Gerreidae | Eugerres lineatus (Humboldt, 1821) | 9 | 6.3-7.5 | 0.44 | 0.64 | 2.67 | 1.25 |
Carangidae | Caranx caninus Günther, 1867 | 10 | 4.0-13.4 | 0.48 | 1.01 | 2.14 | 1.21 |
Centropomidae | Centropomus nigrescens Günther, 1864 | 15 | 4.7-9.6 | 0.73 | 0.64 | 2.14 | 1.17 |
Lutjanidae | Lutjanus argentiventris (Peters, 1869) | 14 | 5.1-10.3 | 0.68 | 1.09 | 1.60 | 1.12 |
Engraulidae | Anchoa lucida (Jordan & Gilbert, 1882) | 15 | 5.3-8.7 | 0.73 | 0.39 | 2.14 | 1.08 |
Lutjanidae | Lutjanus novemfasciatus Gill, 1862 | 10 | 4.2-9.4 | 0.48 | 0.51 | 2.14 | 1.04 |
Mugilidae | Mugil cephalus Linnaeus, 1758 | 2 | 13.0-23.4 | 0.10 | 1.45 | 1.07 | 0.87 |
Ariidae | Ariopsis guatemalensis (Günther, 1864) | 15 | 3.4-8.1 | 0.73 | 0.21 | 1.60 | 0.85 |
Hemiramphidae | Hyporhamphus naos Banford & Collette, 2001 | 5 | 17.4-22.5 | 0.24 | 1.19 | 1.07 | 0.83 |
Poeciliidae | Poeciliopsis fasciata (Meek, 1904) | 14 | 1.6-2.2 | 0.68 | 0.02 | 1.60 | 0.77 |
Engraulidae | Anchovia macrolepidota (Kner, 1863) | 15 | 6.5-14.3 | 0.73 | 1.01 | 0.53 | 0.76 |
Paralichthyidae | Citharichthys gilberti Jenkins & Evermann, 1889 | 7 | 5.8-14.1 | 0.34 | 0.72 | 1.07 | 0.71 |
Clupleidae | Lile nigrofasciata Castro-Aguirre, Ruiz-Campos & Balart, 2002 | 12 | 6.2-7.4 | 0.58 | 0.32 | 1.07 | 0.66 |
Carangidae | Oligoplites saurus (Bloch & Schneider, 1801) | 4 | 10.3-14.5 | 0.19 | 0.62 | 1.07 | 0.63 |
Ariidae | Cathorops liropus (Bristol, 1897) | 3 | 12.2-13.7 | 0.15 | 0.63 | 1.07 | 0.61 |
Lutjanidae | Lutjanus colorado Jordan & Gilbert, 1882 | 2 | 8.9-15.8 | 0.10 | 0.60 | 1.07 | 0.59 |
Centropomidae | Centropomus viridis Lockington, 1877 | 5 | 5.6-11.0 | 0.24 | 0.31 | 1.07 | 0.54 |
Tetraodontidae | Sphoeroides annulatus (Jenyns, 1842) | 3 | 4.5-5.2 | 0.15 | 0.14 | 1.07 | 0.45 |
Ariidae | Cathorops steindachneri (Gilbert & Starks, 1904) | 4 | 4.7-5.5 | 0.19 | 0.06 | 1.07 | 0.44 |
Achiridae | Achirus mazatlanus (Steindachner, 1869) | 3 | 3.1-3.8 | 0.15 | 0.02 | 1.07 | 0.41 |
Poeciliidae | Poeciliopsis pleurospilus (Günther, 1866) | 9 | 1.8-3.2 | 0.44 | 0.05 | 0.53 | 0.34 |
Anablepidae | Anableps dowei Gill, 1861 | 3 | 6.6-11.5 | 0.15 | 0.26 | 0.53 | 0.31 |
Mugilidae | Mugil hospes Jordan & Culver, 1865 | 2 | 8.3-10.2 | 0.10 | 0.23 | 0.53 | 0.29 |
Haemulidae | Haemulopsis axillaris (Steindachner, 1869) | 1 | 10.8 | 0.05 | 0.21 | 0.53 | 0.26 |
Haemulidae | Orthopristis chalceus (Günther, 1864) | 1 | 10.1 | 0.05 | 0.19 | 0.53 | 0.26 |
Labridae | Halichoeres dispilus (Günther, 1864) | 1 | 8.4 | 0.05 | 0.07 | 0.53 | 0.22 |
Gobiidae | Aboma etheostoma Jordan & Starks, 1895 | 2 | 1.7-1.8 | 0.10 | 0.00 | 0.53 | 0.21 |
Characidae | Astyanax aeneus (Günther, 1860) | 1 | 4.6 | 0.05 | 0.02 | 0.53 | 0.20 |
Relative values based on the total number of specimens (2,064), total biomass (17,027.83 g) and the number of sampled sites (13). n: numerical abundance; SL: size ranges in standard length; %n: relative abundance; %B: relative biomass; %F: relative frequency; %IV: relative importance value.
Although the LJB system showed a significant environmental variation between seasons, the fish fauna exhibited their main variations across habitat types. The observed number of species per habitat was 35 in the lagoon, 20 in the inlet and 17 in the channel. Significant differences were found between habitat with the standardized richness [E (S1251)], ranging from 13 taxa (channel) to 26 taxa (lagoon) (H = 6.98, P < 0.05) (Fig. 3a). Differences were found between seasons with the E (S710), ranging from 20 taxa (pre-wet) to 33 taxa (pre-dry), although these were not significant (F = 1.7, P > 0.05) (Fig. 3b).

Figure 3 Rarefaction curves of the expected number of species as a function of sample size calculated by a) habitat and b) season; La Joya-Buenavista lagoon-estuarine system, Mexico. The bars indicate the subsample limits for the number of species and specimens. For the standardized richness values sharing the same letter were not statistically different at P < 0.05 with the non-parametric ANOVA and the Mann-Whitney post-hoc test.
Spatial and temporal variation of fish assemblages
The PERMANOVA results conducted showed that fish assemblages differed significantly at the spatial level (habitat, sampling site; P < 0.01). At the temporal level, differences between seasons were also detected, but not by month or their interactions (Table 3).
Table 3 Results of PERMANOVA based on Bray-Curtis dissimilarities on fish abundance data in response to season, habitat, months, sampling sites, and their interactions in the La Joya-Buenavista lagoon-estuarine (LJB) system (Mexico).
Source of variation | d.f. | MS | F | P |
---|---|---|---|---|
Habitat (Ha) | 2 | 25960 | 11.24 | 0.0001 |
Season (Se) | 3 | 5438.5 | 2.35 | 0.002 |
Sampling site (Si, nested in Ha) | 10 | 3748.1 | 1.62 | 0.008 |
Month (Mo, nested in Se) | 2 | 2397 | 1.04 | 0.426 |
Ha × Se | 6 | 5098.3 | 2.21 | 0.0004 |
Ha × Mo | 4 | 2314.3 | 1.00 | 0.474 |
Se × Si | 30 | 2459.6 | 1.06 | 0.348 |
Residuals | 20 | 2310 |
d.f.: degrees of freedom, MS: mean square; F: pseudo-F statistic value, P: P-values by permutation. Permutated residuals under a reduced model (Type III). The maximum number of permutations = 9,999. Significant P-values (<0.05) are in bold.
According to the SIMPER analysis, the average dissimilarity among habitats ranged from 88.3 to 97.3%. Most of the dissimilarity between lagoon-inlet and inlet-channel was due to the species E. currani (14.7 and 16%, respectively), while Centropomus robalito presented a high contribution (12.8%) in the dissimilarity between the lagoon and channel (Table 4). Regarding the functional guilds in the mouth and lagoon, euryhaline species (
Table 4 Discrimination of 10 fish species based on its relative contribution among habitats by SIMPER analysis; La Joya-Buenavista lagoon-estuarine (LJB) system (Mexico).
Species | Ecological category | Functional group | Spatial distribution | Contribution (%) | ||||
---|---|---|---|---|---|---|---|---|
I | C | L | L vs I | L vs C | I vs C | |||
Eucinostomus currani | EU | OV | 14.66 | 3.52 | 16.01 | |||
Centropomus robalito | EU | PV | 2.35 | 12.79 | 10.02 | |||
Diapterus brevirostris | EU | OV | 7.09 | 6.72 | 5.68 | |||
Lile gracilis | E | ZP | 5.72 | 8.62 | 2.49 | |||
Astatheros macracanthus | SF | OV | 4.88 | 7.62 | 4.43 | |||
Dormitator latifrons | E | DV | - | 8.86 | 6.47 | |||
Atherinella guatemalensis | E | ZB | 7.55 | - | 6.90 | |||
Lutjanus argentiventris | ST | PV | 5.17 | - | 4.78 | |||
Gerres simillimus | EU | OV | 4.65 | 5.86 | - | |||
Mugil curema | EU | DV | 4.58 | 2.93 | 2.99 |
Ecological category: E: estuarine resident, EU: marine euryhaline, ST: marine stenohaline, SF: secondary freshwater. Functional group: DV: detritivore, OV: omnivore, PV: piscivore, ZB: zoobenthivore, ZP: zooplanktivore, I: inlet, C: channel, L: lagoon.
The mean fish size differed among the three habitats (H = 62.26, P < 0.01). Standard length data showed that some of the dominant taxa (A. macracanthus, M. curema, E. currani, D. latifrons) were smaller in the lagoon compared with other habitats according to the Mann-Whitney test. In the case of Lutjanus argentiventris, which was not captured in the lagoon, the smaller specimens corresponded to the inlet. On the other hand, L. gracilis and Atherinella guatemalensis were species whose size does not differ between habitats (Fig. 4).
Relationship between fish assemblages and environmental parameters
The DistLM analysis identified four environmental variables strongly correlated with the fish assemblages when tested in isolation (marginal tests; Table 5). However, in the best solution of the full model (sequential tests, adjusted R2 = 0.48), only habitat structure, salinity and depth together were enough to explain significantly 43% of the total variance. Selected variables by the step-wise procedure were displayed as vectors in the dbRDA ordination plot (Fig. 5). Fish assemblages modeled by four predictor variables revealed two gradients. The first gradient was related to habitat structure and salinity, which kept inlet sites separated from lagoon and channel sites. The second gradient was mainly driven by dissolved oxygen and depth, discriminating against the lagoon sites from the channel sites.

Figure 5 Distance-based redundancy analysis biplot based on linear modelling for fish abundance data, La Joya-Buenavista lagoon-estuarine system, Mexico. Vectors indicate correlations between the four environmental parameters of most significance in the model and the axes. The length and direction of the vectors represent the strength of the relationship concerning a unit circle. Symbols represent an individual sampling event coded by habitat. Black triangles: inlet, gray squares: channel, white circles: lagoon.
Table 5 Marginal and sequential tests result for the distance-based linear models based on fish abundance data in the La Joya-Buenavista lagoon-estuarine (LJB) system (Mexico).
Variable | SS trace | F | P | |
---|---|---|---|---|
Marginal test | ||||
Habitat structure | 11536 | 4.25 | 0.001 | |
Depth | 6417 | 2.17 | 0.036 | |
Dissolved oxygen | 6163 | 2.02 | 0.042 | |
Salinity | 6088 | 1.99 | 0.046 | |
Water clarity | 4375 | 1.38 | 0.209 | |
pH | 3694 | 1.15 | 0.300 | |
Water temperature | 1799 | 0.54 | 0.847 | |
Sequential test | ||||
Habitat structure | 11536 | 4.25 | 0.001 | |
Depth | 6713 | 2.74 | 0.009 | |
Salinity | 4091 | 1.76 | 0.050 |
SS: sum of squares; F: pseudo-F statistic value; P: P-value of permutation test for the relationship between the fish assemblages and environmental parameters. Significant P-values (<0.05) are in bold.
DISCUSSION
The specific richness in the LJB system (48 species) was higher than that described in other ecological studies carried out in the Gulf of Tehuantepec. In the southern Gulf, richness varied between 31 species in the Chantuto-Panzacola system (Díaz-Ruiz et al., 2004) and 40 in Carretas-Pereyra (Velázquez-Velázquez et al., 2008). However, in the first case, comparisons between the fish assemblages should be taken with caution, due to the difference in the sampling effort and method. In the second case, the sampling was more comparable with this study, and differences in richness and composition may be due to less marine influence due to mouth dynamics (Gómez-González et al., 2012). Although it has been pointed out that the use of a single fishing gear results in a biased sample of the assemblage (Clement et al., 2014), in estuarine wetlands the use of cast nets is highly effective for achieving a broad taxonomic representation in several shallow habitats, and with accurate estimate of abundances (Sheaves et al., 2007; Sheaves & Johnston, 2009). However, highly complex habitats should be explored to avoid underrepresenting the cryptic fish fauna.
According to our hypothesis, we found that the fish fauna varied significantly within each geomorphic habitat. When abundances of 48 species at 72 sampling events are ordered on the two primary axes of dbRDA (Fig. 5), a spatial pattern response in fish assemblages on the lagoon, inlet and the channel is identified (Table 4; Fig. 3, 5). The lagoon exhibits a greater richness (35 spp.) with smaller sizes, compared to the inlet (20 spp.) and the channel (17 spp.), suggesting differentiated use of the habitat or the existence of assemblages with their ecological attributes. The above could also be reflected in the percentages of the functional guilds in each habitat (e.g., a greater number of omnivores in the lagoon; Table 4).
In the LJB system, the connectivity among the three geomorphic habitats may be a key attribute that defines the small-scale variation of fish assemblages. Similar spatial patterns were described by França et al. (2009) and Neves et al. (2013) in temperate and tropical coastal lagoons. If the habitats offer differentiated resources in food and shelter, the sequential use of these will favor an efficient use of the lagoon-estuarine system as a whole (Abrantes et al., 2015; Loureiro et al., 2016). Also, the results obtained here are consistent with the theory that the seascape/landscape structure plays a determinant role in the structure of fish assemblages (Kimirei et al., 2013).
Although differences were found at the temporal level, the spatial variations among geomorphic habitats were more significant (Table 3), probably due to a higher number of dominant species residing into the system during the entire annual cycle, and whose movements depend on the availability of resources (Lukey et al., 2006; Herzka et al., 2009). Another explanatory factor could be the relative environmental stability of the system. The significant spatio-temporal interactions detected, habitat × season and between sampling sites, can be indicative of seasonal pulses that define variations in physical and chemical properties within the same ecological or spatial unit (Yáñez-Arancibia et al., 1983; Roselli et al., 2013). However, the lowest relationship found at a temporal scale could be because only one annual cycle was evaluated. Therefore, a study in time series would be required to corroborate these results.
The presence of large areas of complex habitats, like mangroves, is a key factor why fish use lagoon estuarine systems as nurseries. Mangroves support more diverse and abundant small-sized fish assemblages, as they favor greater protection from predators and increase food availability (Nagelkerken et al., 2002; Hindell & Jenkins, 2004; Piko & Szedlmayer, 2007; Green et al., 2012; Neves et al., 2013). The sheltered lagoon and the riverine mangroves that border the channel are suitable habitats for the fish fauna in the LJB system. However, during the dry season, the mangrove root system is not available for fishes; thus, the lack of physical structure may result in a reduction of shelter. Therefore, the highest abundance and richness of species in the lagoon, compared with the other habitats, suggests that this area has a significant ecological value. In this order, Franco et al. (2006) and Verdiell-Cubedo et al. (2013) have reported high abundance values in shallow habitats with small patches of macrophytes or perimetral lagoon marshlands, suggesting that juvenile fish use open-lagoon areas for feeding (Arceo-Carranza & Chiappa-Carrara, 2013; Kathoon et al., 2014). The high abundance of species with detritivore, zoobenthivore and zooplanktivore feeding habits (e.g., Mugil curema, Lile gracilis, Atherinella guatemalensis) were recorded in the lagoon area (Table 4), suggesting that both the physical and biotic attributes of this habitat influenced the population traits of these species.
The proximity of mangroves and mixed substrates also contributes to increasing the complexity of the channel; however, this habitat presented less abundance and richness concerning the inlet, considered of low complexity (Lacerda et al., 2014). However, it has been reported that in tropical estuaries, there are lower values of abundance and richness associated with the mouth area (Neves et al., 2013). Even though, in our study, the inlet is the richest species area of all estuarine channel, which may be due to the occurrence of rare species, mainly of stenohaline affinity (Peralta-Meixueiro & Vega-Cendejas, 2011), during the pre-dry and dry seasons. During the rainy season, freshwater inputs modify the hydrology (salinity, dissolved oxygen, temperature), with a higher number of estuarine species attracted to the beach by a high density of such prey (e.g., Litopenaeus spp.). On the other hand, during the dry season, the beach environment becomes similar to the marine zone, with a lower abundance of fish but a higher richness and dominance of some families (e.g., Carangidae and Sciaenidae).
A meaningful relationship has been founded for several environments between the habitat complexity and the structure of animal assemblages. The ecological paradigm holds that the most complex and physically heterogeneous habitats will support more diverse communities (Stewart et al., 2000). An increase in available refuges due to microtopographic variations of the substrate or physical structures also has been shown to reduce inter and intra-specific competition for space, which potentially increases the richness and abundance of fish (Almany, 2004; Consoli et al., 2018).
Some of the dominant species found in the lagoon and the channel are secondary freshwater and estuarine residents (e.g., Dormitator latifrons and Astatheros macracanthus). The distribution and population structure of this group are determined by a temporal pattern affecting hydrological conditions (Díaz-Ruiz et al., 2004; Velázquez-Velázquez et al., 2008). During the rains, larger organisms were found in the channel, while juveniles of these species occurred in the lagoon during the pre-dry and dry seasons, suggesting a local migration from the streams to mesohaline environments.
In the LJB system, families such as Gerreidae and Mugilidae present greater abundances of larger organisms in the inlet than in the lagoon, suggesting that the sandy bars play an important and specific role for providing habitat. For example, in northern Australia, Blaber et al. (1995) found that the tidal environment, located between the mouth of an open estuary and the adjacent lower marine zone, sustains a relatively complex fish assemblage, composed by dependent juveniles as well as groups of opportunistic species that are move through three habitats (estuary, inshore and offshore waters). The structure of a similar assemblage occurs in the inlet of the LJB system, where juveniles of some species (e.g., Lutjanus argentiventris) were found exclusively in this area, and others of distribution extended by all habitats without variations in their size (e.g., L. gracilis).
Water depth could be a parameter strongly correlated to abundance and size structure in the biotic assemblages of an estuarine system. A general pattern observed shows that larger fishes will be more abundant in deeper habitats, whereas smaller ones will be more abundant in shallow environments (Hutchinson et al., 2014; Mustamäki et al., 2015; Becker et al., 2017). In the LJB system, the greatest depths were recorded in the channel, which may explain the occurrence of larger specimens from families as Lutjanidae, Centropomidae and Cichlidae (Fig. 4); this habitat could be used as a transit or feeding zone in mangroves and marginal lowlands. Hypoxia and low-depth are significant barriers to the presence of larger predatory fish in the lagoon. They affect their ability to move and evade predation, due to physiological stress (Patterson & Whitfield, 2000; Shoji et al., 2005), so greater impact towards recruits will be inflicted by small and juvenile piscivores (Baker & Sheaves, 2009); further evidence for a specific trophic structure for each habitat in the environment (Tecchio et al., 2015).
Although there were seasonal differences, the evidence of fish-habitat associations for many species was broad regardless of seasonality, probably because a high percentage of resident species in the system was found during the annual cycle. However, sampling needs to be expanded to compare contiguous systems and between-year variations in order to corroborate these conclusions. Several anthropogenic pressures are affecting the lagoon estuaries of the Gulf of Tehuantepec, whose vulnerability has increased in recent years as a result of the lack of planning in the coastal basin. Due to their ability to respond to environmental changes, fish assemblages are a potential tool for monitoring the ecological integrity of coastal systems, so their study should be continued. Likewise, coastal management should be developed to conserve the mosaic of habitats that integrate all the system, rather than distinct habitats.