Complementary Sampling Methods to Improve the Monitoring of Coastal Lagoons

: Monitoring the ecological status of marine coastal lagoons requires the integration of multiple indices. However, the efficacy of monitoring programs is complicated by the diverse array of habitats that conform coastal lagoons. In this study, we compared four sampling methods (25-m and 50-m beach seines, beam trawl and Riley push net) in the Ria Formosa coastal lagoon (South Portugal) for assessing fish assemblage and diversity. We compared species richness and assemblage structure with species accumulation curves and multivariate analysis, and assessed diversity patterns using taxonomic, phylogenetic and functional diversity indices. Variation in fish assemblage structure was mostly explained by gear type, and almost all diversity metrics varied not only according to sampling method but also depending on habitat characteristics and season. Based on operational costs and diversity patterns captured by each gear, we conclude that the combined use of 25-m beach seine and beam trawl is the preferred approach. The proposed methodology will provide the data necessary for assessment of ecological status of coastal lagoons.


Introduction
Estuaries and sheltered lagoons contain some of the most productive coastal marine habitats such as seagrass beds and salt marshes [1,2]. These ecosystems are essential feeding and nursery grounds for juvenile fishes and invertebrates, including many species with commercial and recreational value [3][4][5][6][7]. However, their long-term monitoring is particularly challenging due to the heterogeneous spatial distribution of species across different habitat types, coupled with seasonal shifts in assemblage structure and diversity [8,9]. Yet, their effective monitoring is fundamental, as an increasing number of anthropogenic threats (e.g., pollution, habitat loss, sea level rise, and overfishing) are degrading the state of these valuable ecosystems [10,11]. Further, in the context of the Water Framework Directive (WFD), member states must establish monitoring programs to provide information on long-term changes for each water body type (lakes, rivers, transitional waters, coastal waters) [12]. Fishes are an important biological component of coastal lagoons, performing key ecological roles that underpin ecosystem productivity and resilience (e.g., top-down control of prey populations, sediment reworking, and nutrient recycling) [13][14][15]. Due to their importance, changes in fish assemblage metrics (e.g., diversity, composition and abundance) are part of the quality elements for the assessment of the ecological status of different water bodies under the WFD.
Until recently, the assessment of fish assemblages in coastal lagoons has mainly focused on traditional taxonomic diversity indices (TD, e.g., species richness, Shannon-Wiener diversity, and Pielou's evenness) [16]. However, these metrics treat all the species equally without considering their potential contribution to a range of ecosystem functions. In contrast, functional (FD) [17] and phylogenetic (PD) [18] diversity metrics incorporate information on species ecological similarities based on traits (i.e., morphological, ecological, physiological, and behavioural characteristics) and evolutionary information [19,20]. There is growing empirical evidence, in both terrestrial and marine environments, of the importance of these often-overlooked biodiversity components for ecosystem functioning and resilience [21][22][23].
Given the heterogeneity of habitats and seasonal variability in coastal lagoons, sampling their fish assemblages with the appropriate gear or combination of gears is essential for obtaining reliable data for the application of fish-based ecological indicators [24][25][26]. There are a variety of fishing gears for sampling coastal lagoons, so understanding their biases and whether they can replace or complement each other is crucial [27,28]. Beam trawls are fishing gears commonly used not only for commercial fishing activities but also for regular sampling and monitoring of fish communities in estuaries [3,24,29]. This active gear is suitable for sampling large numbers of fish, especially demersal species [30]. Other methods such as beach seines and push nets require different physical characteristics of the study site, since both need access from the shore to shallow waters. Trawls and seines frequently have variable catch efficiency, due to either differences in gear design or gear avoidance by certain fish groups [31]. This sampling bias, together with the fact that estuaries contain highly variable hydrographic and spatial-temporal characteristics, justify the need of multi-method surveys to capture an accurate representation of fish assemblage structure and diversity. Yet, few studies have examined to what extent the inherent bias associated with different sampling methods might undermine the ability of monitoring programs to detect shifts in important but often overlooked functional and phylogenetic diversity components [32].
In the present study, we compared four commonly used sampling methods for coastal lagoon systems-beach seines of 25-m and 50-m, beam trawl and push net-to assess their complementarity in fish assemblage metrics (composition and structure, taxonomic diversity, functional diversity and phylogenetic diversity). Specifically, we hypothesize that these four gears capture different representations of species assemblage and diversity patterns. Other studies have compared different fishing gears to determine their efficiency in sampling species composition and abundance [27][28][29]33], but this study is the first to assess differences among this particular combination of gears using less conventional phylogenetic and functional diversity indices. The results of this work provide useful information for management agencies aiming at identifying which method or combination of methods is more suited to track changes in the ecological status of these valuable but highly threatened ecosystems.

Study Site and Sampling Stations
The Ria Formosa is a tidal lagoon extending about 55 km along the south coast of Portugal, consisting of salt marshes, subtidal channels and tidal flats covering a surface area of approximately 170 km 2 , up to 6 m in depth [4,9]. With tidal elevations of 1.30 and 2.80 m at mean neap tide and spring tide, respectively, the minimum and maximum areas covered by water during spring tides are 14.1 and 63.1 km 2 [34]. Since it is located between the sea and the land, this lagoon has distinctive biological, ecological and hydrodynamic features, with a variety of different habitats that can be distinguished based on substrate type (sand, gravel or fine mud), depth, vegetation and distance to the sea [4]. The lagoon has extensive patches of seagrass (Cymodocea nodosa, Zostera marina and Zostera noltii) providing shelter from predators and adverse weather conditions, as well as food sources [9].
Approximately 90% of the area is a natural park corresponding to category V in the International Union for Conservation of Nature (IUCN) classification of protected areas [35,36] and has high ichthyological diversity compared with other equivalent ecosystems [37,38]. This area was integrated as a Natura 2000 site under the Birds and Habitats European Directives [39] and it was also declared a RAMSAR site for conservation of wetlands in 1981 [40]. Nevertheless, several economic activities take place (e.g., fishing, harvesting of bait, aquaculture, tourism, shipping, airport activity), putting pressure on this vulnerable ecosystem [41].
In order to sample the fish fauna in the different habitats present in the Ria Formosa, a variety of sampling methods were tested, and the choice of sampling stations was based on a stratified sampling strategy. Sampling stations were chosen by first dividing the lagoon into three areas with very distinct hydrodynamic characteristics: (1) areas of strong coastal influence, near the inlets and with a strong hydrodynamic regime influenced mostly by coastal waters (closer to barrier islands); (2) interior areas with salt marshes, which correspond to shallow vegetated areas greatly affected by the terrestrial environment and usually dry at low tide; and (3) main and secondary channels, that represent the deepest parts of the Ria Formosa. For these three main areas, 59 sampling stations were chosen based on the types of habitats present and the sampling gear to be used ( Figure 1). Habitats were classified according to the presence or absence of vegetation (VEG/UNVEG). Seasons were defined as autumn (AU, from October to December), winter (W, January to March), spring (SP, April to June) and summer (S, July to August). Monthly samples were collected from September 2000 to April 2002, with 37 stations chosen for the 25-m beach seine, 4 stations for the 50-m beach seine, 12 stations used for the beam trawl, and 6 stations for the Riley push net ( Table 1). The selection of sampling stations was based on visual surveys of the area and existing aerial photographs and maps [4]. Due to logistic problems, 92% of the target sampling was achieved for the beam trawl, Riley push net and 50-m beach seine, and 84% for the 25-m beach seine. For the latter sampling gear, the main factors affecting sampling were the short period for sampling (could only be done at low tide), the large number of locations, the distances between them, and occasionally, problems such as the net getting stuck, requiring repeating the sample.  The two beach seine nets were deployed at low tides in the margins of the main channels only when the amplitude of the tide was less than 2 m, and during a period of 2 h before and 2 h after the low tide to reduce variability in tidal amplitude between sampling events. One of the beach seine nets was 25-m long, 3.5 high in the middle and was made of 9-mm mesh netting. The net was towed by a boat from one end and by two people on shore from the other end, resulting in an average sampled area of 1087 m 2 (based on GPS measurements). The 50-m beach seine was 3.5 m high in the middle and was made of 13 mm mesh [42]. For the setting of the net, one end was held on shore while the net was set in a circle; no towing took place. Based on GPS measurements, the average area sampled was 295 m 2 , but since 3 replicates were taken at each station and the samples pooled across replicates, the average total area sampled was 885 m 2 for each station.

Beam Trawl
A beam trawl 2.6 m wide and 0.45 m high at the mouth was used at low tides in the main channels. The cod end was 10 m long and made of 9-mm mesh netting. Tows of 300 m were performed at 1 knot, resulting in a swept area per tow of 780 m 2 .

Riley Push Net
The Riley push net was 1.5 m wide and 0.5 m high at the mouth, with a cod end split into two "trouser" legs of 2-mm fine mesh netting. Three 30-m samples were taken at each sampling station within interior areas of the Ria by one person who stood between the two legs of the cod end, resulting in a total swept area of 135 m 2 .

Processing of Fish Collected
The same method of processing samples was applied for all gears. Invertebrates and species such as seahorses (Hippocampus hippocampus and H. guttulatus) were caught and released alive, while the other fish species were placed immediately after capture in boxes with an ice flurry to minimize suffering and transported to the laboratory for processing. At the laboratory, the fish were sorted, identified to the lowest taxonomic level possible and counted. Total length was measured to the nearest mm. Detailed information about the sampling gears and habitats sampled together with photographs are available in Erzini et al., (2002) [4] (report available upon request).

Assemblage Composition and Structure
A Euler diagram was constructed in order to compare the overlap in species composition caught by beach seine (25-m and 50-m included in the same group), beam trawl and Riley push net. To estimate the relationship between the number of species observed and the sampling effort as a measure of sampling efficiency of the different gears used, species accumulation curves were computed via sample-based rarefaction [43]. The method calculates the expected species richness for each sample under random order from 1000 permutations of the data, and sampled area was used to standardize effort between sampling methods. For very high sampling efforts, the curves would eventually reach an asymptote matching the assemblage richness available to the method chosen and the more concave-downward the curve, the better sampled the community [44].
A multivariate regression tree (MRT) was constructed to explore the relative influence of sampling gear, season and habitat in explaining fish assemblage structure. MRT is a statistical technique that combines multivariate regression and constrained clustering, since it forms clusters based on a measure of species dissimilarity that are defined by a set of predictor variables [45]. In this case, fish abundance data were partitioned successively in two subsets by selecting one of the factors (gear, season or habitat) that maximizes the homogeneity of the resulting clusters [46]. Each cluster defines a species assemblage and is determined by the associated explanatory factors (gear, season or habitat); this procedure is graphically represented by a tree with nodes where the groups are split and terminal nodes define the final clusters. Optimal tree size was chosen by minimization of cross-validated relative error (CVRE) and smallest tree size. MRT analysis was chosen due to its ability to deal with categorical variables and highorder interactions among explanatory variables [45]. Discriminant species were identified within the tree as species that contribute most for the explained variance at each node. Indicator values (IndVal) were calculated for each species in each terminal group by multiplying a measure of specificity ( , mean abundance of species i in the sites of group j compared to all groups in the study) and a measure of fidelity ( , relative frequency of occurrence of species i in the sites of group j) [47]. This index ranges from 0 (no occurrences of the species within a group) to 1 (the species occurs at all sites within the group and does not occur at any other site). Species with high index values (≥ 0.2) for a cluster were considered indicator species [48]. Species abundance was converted to densities (numbers per sampled area) and standardized by dividing species density at each sampling station by the total density for all species at that same station.

Diversity Indices
For the traditional taxonomic diversity analysis, the following metrics were computed for each sample of each gear: Shannon-Wiener diversity index (H), Pielou's measure of evenness (J) and species richness (S). Species richness represents counts of the number of species (S). The Shannon-Wiener index is defined as . Abundance data was converted to densities (number per sampled area) before calculation of the indices and square root transformation used to balance the contribution of dominant and rare species.
As a metric of phylogenetic diversity, the index of taxonomic distinctness (∆*) was chosen and is defined as the average taxonomic path length between two individuals chosen at random from the sample, traced through a standard Linnean classification tree, conditional on them being from different species [49]. It is calculated by dividing taxonomic diversity (∆) by the Simpson index and is not affected by the evenness properties of the species abundance matrix. This phylogenetic diversity index takes the form: ∆ * = ΣΣ (ΣΣ ) , where ω represents taxonomic distances among taxa, x are species abundances and the double summation goes over species i and j [50]. Taxonomic distances were calculated among taxa at variable step lengths based on an aggregation matrix by species, genus, family and order. Species abundance matrices were transformed to square root density data and ∆* calculated for each sample. The functional diversity (FD) was measured by the Rao index of functional diversity which represents the probability that, choosing two random species from the sample, they have different trait values or trait categories [51]. The Rao index combines matrixes of species abundances with matrixes of species dissimilarities based on traits' differences among species [52] and is defined as = ∑ ∑ , where S is the number of species, pi and pj are the proportion of ith and jth species and dij the dissimilarity between species i and j. Following Bosch et al., (2017,2021) [32,53], maximum body length, trophic breadth and trophic level were included as continuous (scaled between 0 and 1), while trophic group, water column position, preferred substrate and body shape were considered as categorical traits, using fuzzy coding for traits with more than 2 categories (Table 2). Information on traits was collected from the published literature and Fishbase [54], and when species specific attributes were not available, values from species within the same genus and geographic range were used. The Rao index was calculated using the Macro excel file "FunctDiv.exl" [51] on density data (numbers per sampled area), first for each trait and then averaged for each sample across all traits together. We tested for differences in diversity indices between sampling gears, season, habitat characteristics (vegetated and unvegetated), and main area using generalized linear models (GLM), applied for each diversity index independently. Second order interactions were included between gear-vegetation, and gear-season. The area consisted of 3 factors: inner, channels and outer areas (I, C and O in Figure 1), and there was no interaction included between gear-area since not all gears were used for all the three areas. The following error distributions were fitted: Poisson distribution for species richness (S) with log link function; normal distribution to species diversity (H), species evenness (J) and functional diversity (FD) (identity link function); and gamma distribution to phylogenetical diversity (PD) (identity link function).
All analyses were carried out using R statistical software version 3.6.0 [56]. The packages 'vegan' [57] and 'biodiversityR' [58] were used for computing species accumulation curves and calculation of diversity indices (except for the Rao index of functional diversity). The multivariate regression tree (MRT) and species indicator values were calculated with the packages 'mvpart' and 'MVPARTwrap' [59]). GLMs were conducted with the 'stats' package from base R for model evaluation and residuals analysis.

Assemblage Structure
During this study, a total of 255,345 fish belonging to 106 species within 33 families (103 teleosts and 3 chondrichthyes) were captured. Atherina spp., Sardina pilchardus, Gobius niger and Pomatoschistus microps represented 70% (range of 42 to 94%) of the total catches across gears (Supplementary Materials, Table S1). In terms of species composition, there was considerable overlap between the different gears (40 species; Figure 2). Beach seines and beam trawl shared the highest number of species (23 species), while beam trawlpush net and beach seines-push net shared only 2 and 3 species, respectively. Beach seines accounted for 29 species that were not caught by the other gears, mainly species of the Sparidae, Labridae and Triglidae families. The beam trawl caught eight species that none of the other gears captured, with Soleidae being the most representative family (>85% abundance). Two species were only caught by Riley push net, Dentex macrophthalmus and Parablennius sanguinolentus, with forty-five species caught in common with the other gears. Overall, only the 25-m BS curve was close to reaching the asymptote, having by far the highest cumulative sampled area ( Figure 3). Nevertheless, for comparable sampling efforts (e.g., less than 0.05 km 2 ), both the 50-m beach seine and beam trawl caught more species than the 25-m beach seine. The species accumulation curve for the Push net is difficult to distinguish from the others due to low sampling effort.
The best multivariate regression tree explaining variation in fish assemblage structure consisted of four nodes and five leaves (terminal nodes; Figure 4). This model explained 27.22% of the variation in species abundances. The first split in the tree explained 15.87% of the variability and separated the fish assemblage according to two combinations of fishing gears-beam trawl-push net, and 25-m and 50-m beach seine, for which the main species contributing to explaining this split belonged to the genus Atherina spp. The second node separated the assemblage between the other pair of sampling gears (push net and beam trawl). For the 25-m BS-50-m BS cluster (right side of node 1), the species assemblage was divided according to season, with Summer (S) containing a distinct assemblage when compared with the rest of the seasons. On the left side of the main node, the last split between Habitat types explained 4.57% of variation, with habitats containing vegetation (VEG) clearly distinct from non-vegetated (UNVEG).  Overall, five distinct assemblages were identified, each one with a different set of indicator species. The first and second clusters were defined by samples obtained only with the Riley push net; the first group (I) is dominated by the common sand goby (Pomatoschistus microps) in sandy/muddy grounds; and the second group (II) contains three species of pipefish (Syngnathus typhle, S. abaster and Nerophis ophidion), Baillon's wrasse (Symphodus bailloni), rock goby (Gobius paganellus) and the two-spotted clingfish (Diplecogaster bimaculata) in seagrass habitats. The third assemblage (III) is defined by species captured with beam trawl and the analysis did not separate at the habitat level, with toadfish (Halobatrachus didactylus), two goby species (black goby, Gobius niger, and sand goby, Pomatoschistus minutus), two seahorse species (Hippocampus guttulatus and H. hippocampus), the grey wrasse (Symphodus cinereus), the small red scorpionfish (Scorpaena notata), and the flatfish Arnoglossus thori representing the indicator species for this assemblage. The last two groups (IV and V) were classified by samples collected both with the 25-m and 50-m beach seines, and the only split was defined by season. The summer assemblage was dominated by five sea bream species (black seabream, Spondyliosoma cantharus, Senegal seabream, Diplodus bellottii, White seabream, D. sargus, two-banded seabream, D. vulgaris, and gilt-head seabream, Sparus aurata), European pilchard (Sardina pilchardus), striped red mullet (Mullus surmuletus) and the European bass (Dicentrarchus labrax). In the other seasons, only the resident silverside species was identified as indicator species (Atherina spp.).

Diversity Indices
There were clear differences in all diversity measures between sampling gears, presence of vegetation and season ( Figure 5, Table S2 in Supplementary Materials). In terms of species richness, the highest values were for the 50-m beach seine (50-m BS) and beam trawl (BT) (panel A). Habitats with vegetation sampled with any of the four gears showed higher species richness than unvegetated locations. There were higher values of species richness (S) for 25-m BS, 50-m BS and beam trawl during summer months, and lower values for all gears in winter (panel B in Figure 5; Table S2). For species diversity represented by the Shannon-Wiener index (H), the highest values were registered for 50m BS and beam trawl (panels B and C). There were also higher values of H for vegetated locations compared to unvegetated locations, and higher values of H for stations sampled during summer months (particularly 50-m BS), and lower values of H during the winter period. There is a distinct effect of sampling vegetated habitats with push net, with high values of H when compared with habitats without vegetation sampled with 25-m BS (pvalue < 0.0001; Table S2). With regard to species evenness (J), there were slightly higher values for BT, and very low values for samples collected with push net in unvegetated locations (panels E and F). There was not an obvious effect of season, only slightly lower values for samples collected in spring with 25-m BS, and samples collected by push net had high variability. In terms of phylogenetic diversity (PD), 25-m BS had highest values (panels G and H). The GLM detected a significant interaction between push net and vegetation, meaning that vegetated habitats sampled with push net had higher phylogenetic diversity than unvegetated stations with 25-m BS. For functional diversity (FD), there were no major differences between sampling gears (panels I and J). There were higher values of FD in vegetated habitats sampled with any of the four gears, and there was a strong effect of sampling with push net in vegetated habitats compared with unvegetated. Summer and spring months show slightly higher values of FD but this is not very clear for samples collected with 25-m BS and 50-m BS.

Discussion
For the development and implementation of monitoring programs in estuaries and coastal lagoons, it is necessary to select the most adequate sampling methods for detecting spatiotemporal changes in species composition, abundance, and diversity, while also minimizing the costs and damage to the local assemblages of these sensitive ecosystems [28,32]. Here, we showed that the use of complementary sampling gears suitable for particular habitat types within a coastal lagoon capture a wide range of ichthyofaunal diversity that has been linked to the ecological status and the functioning of these highly productive ecosystems [60].
A strong seasonal influence was found in previous studies in the Ria Formosa, where the fish assemblage richness and abundance increased with the recruitment of marine juvenile migrant individuals during spring and summer into the lagoon [42,61,62]. This might explain the identification of two assemblages distinguished by season in the MRT and a significant effect of season in the GLMs in terms of species richness and diversity. Furthermore, different habitats had distinct characteristics not only in terms of type of substrate, but also depth at low tide, hydrodynamics and distance to the openings of the sea, environmental factors that play a strong role in structuring ichthyofauna diversity and abundance in the Ria Formosa [4]. For example, Ribeiro et al., (2012) [63] found that sampling with Riley push net in vegetated habitats resulted in higher species richness and diversity than unvegetated habitats, where higher densities of fewer species were encountered. This gives support to the results, where a distinct assemblage was identified (cluster II in MRT and significant interactions in the GLMs) characterized by species living in vegetated habitats that were sampled with the push net. Similar outcomes of seagrass and saltmarsh habitats containing maximum diversity and richness were found in other studies [9,37].
Not all diversity metrics (species richness, species evenness, Shannon entropy, taxonomic distinctness, and Rao's quadratic entropy for functional diversity) varied significantly between sampling gears (Table S2). The results show no significant differences in terms of functional diversity between 25-m beach seine and beam trawl, despite the latter having higher values of species richness and diversity (p-value > 0.05). This could be because of a large range of FD values across samples, especially for the 25-m beach seine. This variation can also be explained by the nature of each gear itself; some gears are more efficient in sampling particular habitats or species with particular functional traits than others. The beach seines are more efficient in sampling small pelagics since they fish the entire water column, particularly the 50-m beach seine (highest species accumulation rate, Figure 3), while the beam trawl catches mostly benthic and epibenthic species since it operates close to the seafloor. The Riley push net captures a large portion of small individuals, being most efficient for sampling juveniles of commercial species [63]. In fact, the beam trawl was used to sample the main deeper channels in the lagoon and captured mainly benthic species, while the beach seines were deployed in the secondary channels and caught more benthopelagic species, belonging to Sparidae and Labridae families. In addition, the deployment methods of the two beach seines differed, with the 50-m used to encircle while the 25-m net was first towed along the shore before being hauled to land. However, the analysis did not show significant differences between these two gears.
Each index gave a distinct picture of the fish assemblage, with only a few similar comparisons between gears across indices. This shows the importance of considering alternative measures of diversity [32]. Other studies have found that functional diversity can be positively related with species richness and diversity [64], but this relationship is not always positive and linear [52]. Phylogenetic diversity and species richness also display different patterns of diversity and do not seem to be correlated [65,66]. Functional traits may change with species ontogenetic shifts, especially in areas that contain different life stages [67]. However, information on trait variability for the life stages of all species included in the analysis was not available but should be accounted for in future studies. Indices based on biomass can be useful complements to those based on abundance, especially for revealing different insights into temporal trends of functional diversity [68], but this is outside the scope of this study. Following the Water Quality European Standard from 2006, the recommended methodology for sampling species composition and abundance in generic transitional waters such as the Ria Formosa is with the beam trawl [69]. However, as found in our comparative study, other sampling gears such as beach seines and push net provide different representations of the fish assemblage, particularly when different dimensions of diversity (FD and PD) are incorporated. Even though the beam trawl samples contained higher species richness and diversity, they had lower phylogenetic diversity. For sampling with a push net, the fish assemblage changed drastically between habitats with and without vegetation and could only be operated at specific locations. As such, a combination of 25-m beach seine and beam trawl might be an effective sampling strategy to cover multiple aspects of diversity. As in Rotherham et al., (2012) [33], each gear used (beam trawl and multi-mesh gillnets) gave a unique picture of assemblages of fauna, so the most complete representation of the entire fish community required sampling with both methods. Several multi-metric fish indices use distinct sampling gears [25]. In Ireland for example, a standard multi-method approach (beach seine, fyke nets and beam trawl) is used in transitional waters for the implementation of the WFD and the development of an estuarine multi-metric fish index [70].
When designing a monitoring plan, the relative costs of deploying each sampling gear need to be taken into account (e.g., human resources, environmental impacts, time and financial expenditures). For this study, all the sampling gears required transport by boat to the sampling locations. The Riley push net required only one person to operate, while the beam trawl needed the skipper to navigate the boat and deploy the gear, and at least another person to help record data, label and store the catches. The beach seines demanded more human power (skipper and three people) to haul the nets and process the catches, particularly the 50-m beach seine. Additionally, the beam trawl operates with heavy ground gear that drags along the lagoon floor and disturbs bottom habitats. On the other hand, the beach seines need to be hauled by a group of people that were occasionally stepping on the seagrass patches which can damage these sensitive habitats. Although operating the beach seines was more labour intensive, these gears allowed sampling of a greater variety of habitats. In contrast, the beam trawl and push net were limited to the deeper channels and the shallow creeks, respectively.

Conclusions
Given the heterogeneity of habitats, variability among sampling gears, and seasonal effects, the use of a multi-gear approach would provide a robust assessment of the fish assemblage structure in coastal lagoons as the Ria Formosa. Combining the 25-m beach seine and beam trawl might be the most advantageous strategy given the limitations of sampling with a push net and the operational costs of the 50-m beach seine. This work thus contributes with new knowledge that adds to current guidance on selection of fish sampling methods in coastal lagoons, an essential parameter for the assessment of ecological status and biodiversity conservation of these ecosystems. This information is of paramount importance for implementing policies and management plans at local, regional, and national level to meet the objectives set out in the UN Sustainable Development Goals under the 2030 agenda.
Supplementary Materials: The following supporting information can be download at https://www.mdpi.com/article/10.3390/d14100849/s1, Table S1: All species identified and counted for each gear; Table S2 Funding: This study received funds from the Commission of the European Communities (DG XIV C1/99/061) and Portuguese national funds from FCT -Foundation for Science and Technology through projects UIDB/04326/2020, UIDP/04326/2020 and LA/P/0101/2020.

Institutional Review Board Statement:
At the time when the fieldwork was conducted (2000)(2001)(2002), there was no code of practice to handle live animals, but the handling practices included the 3R approach ("Replacement, Reduction and Refinement") through subsampling highly abundant fish species, releasing of animals in good condition and practices for minimizing suffering and pain of animals. Permits authorizing the sampling of fish in the Ria Formosa lagoon were obtained from the Institute for the Conservation of Nature and Forests (ICNF), the Ria Formosa Natural Park (PNRF) and Port of Faro Maritime Authority (CPF/CPO/IPTM).

Data Availability Statement:
Data requests should be addressed to the corresponding author.