Spatial Patterns of Macrozoobenthos Assemblages in a Sentinel Coastal Lagoon: Biodiversity and Environmental Drivers

: This study presents an assessment of the diversity and spatial distribution of benthic macrofauna communities along the Moulay Bousselham lagoon and discusses the environmental factors contributing to observed patterns. In the autumn of 2018, 68 stations were sampled with three replicates per station in subtidal and intertidal areas. Environmental conditions showed that the range of water temperature was from 25.0 ◦ C to 12.3 ◦ C, the salinity varied between 38.7 and 3.7, while the average of pH values ﬂuctuated between 7.3 and 8.0. In vegetated habitats, biomass values of the seagrass Zostera noltei Hornemann ranged between 31.7 gDW/m 2 and 170.2 gDW/m 2 while the biomass of the seagrass Ruppia cirrhosa (Petagna) Grande between 54.2 gDW/m 2 and 84.7 gDW/m 2 . Sediment analyses showed that the lagoon is mainly composed of sandy and silty sediments. We recorded 37,165 individuals of macrofauna distributed in 63 taxa belonging to 50 families, with a mean abundance value of 4582.8 ind/m 2 and biomass average of 22.2 g/m 2 . Distance-based linear modeling analysis (DISTLM) identiﬁed sediment characteristics, water parameters and habitat type (biomass of Z. noltei ) as the major environmental drivers inﬂuencing macrozoobenthos patterns. Our results clearly revealed that the hydrographic regime (marine and terrestrial freshwater), sediment distribution and characteristics and the type of habitat (vegetated vs. unvegetated substrate) are the key factors determining the species composition and patterns of macrozoobenthos assemblages.


Introduction
Coastal lagoons are among the marine habitats with the highest biological productivity [1] and perform an important ecological function by providing forty-one varieties of goods and services [2]. However, coastal lagoons are semi-enclosed coastal systems (SECS) where environmental conditions are highly changeable due to their confined nature and their shallowness. SECS are especially vulnerable to the impacts of human activities resulting from mining, industry, tourism and urban development [3,4]. The geomorphology of these SECS renders them particularly vulnerable to global changes, such as sea-level rises, increased temperatures, storminess, droughts, floods and changes in sediment dynamics. They are "hotspots" of global change and vulnerability to environmental, economic and social pressures [5].
Coastal lagoons are sentinel systems that are highly vulnerable to the impacts of climate change [6]. They have natural conditions that play a key role in regulating water structure and spatial patterns. This study also aims to highlight the environmental drivers that govern the spatial distribution of benthic communities.

Study Area
The Moulay Bousselham lagoon is the northernmost lagoon on the Moroccan Atlantic coast (Figure 1). The lagoon is located 125 km north of Rabat; it has an elliptical shape, with a maximum length of 9 km, a maximum width of 5 km and an area of 35 km². The communication of the lagoon with the Atlantic Ocean is done through a narrow, sinuous and relatively deep gully (up to 6 m), which branches out in the direction of the lagoon by shallow subtidal channels, ensuring the circulation of water during the flood and ebb. The freshwater supply is provided by two rivers: Canal Nador in the south and Oued Drader in the northeast of the lagoon. The tidal part of the latter course divides the lagoon into two body-waters locally known as 'Merja': (i) The Merja Kahla, which is extended on 3 km on the north part of the lagoon, is very shallow, and its bottom is covered with a very dark mud; and (ii) the Merja Zerga, which represents the major part of the lagoon as it is extended over 27 km² and appears blue due to its high depth at high tides [21]. The depth of the lagoon varies between 0 and 2 m depending on the tidal cycle and rainfall. During the annual cycle, the average salinity of the lagoon water fluctuates from 24.0 to 36.3 at high tide and from 8.0 to 32.5 at low tide [25].

Sample Collection and Environmental Analyses
Grid sampling design encompasses the entire intertidal and subtidal areas of the Moulay Bousselham lagoon, with a combination of sample points taken at 500 m intervals. Sixty-eight stations ( Figure 1) were sampled in the autumn of 2018 with three replicas per station. In subtidal areas, the samples were collected using a Van Veen grab, and each sample had a surface area of 0.1 m². While in intertidal zone, samples were taken using a PVC corer with a diameter of 12.5 cm, and each replica was a fusion of 10 cores, covering a total area of approximately 0.12 m².

Sample Collection and Environmental Analyses
Grid sampling design encompasses the entire intertidal and subtidal areas of the Moulay Bousselham lagoon, with a combination of sample points taken at 500 m intervals. Sixty-eight stations ( Figure 1) were sampled in the autumn of 2018 with three replicas per station. In subtidal areas, the samples were collected using a Van Veen grab, and each sample had a surface area of 0.1 m 2 . While in intertidal zone, samples were taken using a PVC corer with a diameter of 12.5 cm, and each replica was a fusion of 10 cores, covering a total area of approximately 0.12 m 2 .
The samples were sieved in situ through a 1 mm mesh. The material retained on the mesh was fixed and preserved in seawater with formalin (4%) and colored with Rose Bengal. In the laboratory, macroinvertebrates were sorted, identified and counted. Biomasses were obtained after calcination in the oven at 450 • C for 4 h.
Physicochemical parameters (water temperature, salinity and pH) were also measured in situ with a HANNA portable multiparameter. Each sample of the macrofauna was accompanied by an additional sediment sample to determine their precise granulometry, carbon content and total organic matter (TOM).
Grain size was measured using a laser particle size analyzer (Malvern Mastersizer 2000©) after preparing the sediments in a sodium hexametaphosphate solution [26]. The grain size distribution was then treated with the Gradistat © Excel package [27]. Mean, sorting, skewness, kurtosis, decile statistics (including the median used to characterize the sediment type: d50) and clay/silt/sand composition were calculated to precise the textural group of each sample [28,29]. A LECO© carbon analyzer estimated the carbon, CO 2 and CaCO 3 percentages after 1400 • C dioxygen burning and mineral decarbonizing with sulfuric acid solution [30], while the organic matter content was determined by estimating the total organic matter (TOM). The samples were oven dried (at 60 • C for 48 h) and ignited in an oven for 4 h at 500 • C. The percentage weight loss during the ignition step is reported as TOM [31]. When present, the biomasses of the seagrasses Zostera noltei Hornemann and Ruppia cirrhosa (Petagna) Grande were measured using a dry weight (gDW/m 2 ). The seagrasses were isolated and rinsed with water and then dried in an oven at 60 • C for 48 h.

Data Analysis
The data matrix with macrofaunal abundance per station was transformed into square root, and then the Bray-Curtis similarity was calculated between stations. The similarity matrix was analyzed using agglomerative hierarchical clustering (AHC) to identify macrofaunal affinities. The environmental variables were transformed to Log(X + 1). Percentage similarity analysis (SIMPER) was used to identify the taxa that contributed most to disparities between each identified assemblage and to the dissimilarity among them.
DISTLM analysis (Distance-based linear modeling) was used to assess the contribution of environmental variables to the variability observed in the macrofaunal assemblages [32]. Results were visualized using the graphical representation of the ordination method of redundancy analysis (RDA). RDA is a constrained ordinate used to identify the linear combinations of predictor variables that explain the greatest variation in the species/abundance matrix, i.e., it shows the pattern of species/abundance (response) data as constrained by the predictor variables [33].
Spatial distribution and biodiversity were described by univariate analyses based on the following parameters: abundance (N, the number of individuals per m 2 ), species richness (S), Shannon-Wiener diversity index (H ) [34] and Pielou's evenness index (J ) [35]. All these analyses were performed with PRIMER v6.0 software [36,37].
The (S) (log 2 S) and (H ) indices of the assemblages were plotted together on a twodimensional graphical representation in the Diversity Model (DIMO) considered as a synthetic tool [38].

Environmental Variables
The spatial variation of environmental conditions is shown in Figure 2. The range of water temperature was from 12.3 • C at station I50 to 25.0 • C at station S1. The salinity varied among stations, with the maximum recorded at station S2 (38.7) and the minimum at the station I50 (3.7). The pH values fluctuated between 7.3 (station I9) and 8.0 (station I22). In stations with vegetated habitats (Figure 2), biomass values of Zostera noltei seagrass ranged between 31.7 gDW/m 2 in station I29V and 170.2 gDW/m 2 in station I23V, while the biomass of Ruppia cirrhosa seagrass ranged between 54.2 gDW/m 2 in I6V and 84.7 gDW/m 2 in I13V.

Benthic Macrofauna
Overall, 37165 individuals of benthic macrofauna were collected and are distributed on 63 taxa including 50 families. Mollusca was the predominant phylum with 24 species belonging to 18 families. For Arthropoda, 20 species were counted with 19 families. The phylum Annelida was the third dominant group with 18 species belonging to 12 families. One family represented the phylum Nemertea.

Benthic Macrofauna
Overall, 37165 individuals of benthic macrofauna were collected and are distributed on 63 taxa including 50 families. Mollusca was the predominant phylum with 24 species belonging to 18 families. For Arthropoda, 20 species were counted with 19 families. The phylum Annelida was the third dominant group with 18 species belonging to 12 families. One family represented the phylum Nemertea.
The DIMO model distinctly separated the community groups and displayed a type 4 dynamic (non-constant type), where all three parameters (S, H and J ) changed ( Figure 6). According to the DIMO model, the stations located in the muddiest areas (near the Nador Canal), some stations in the subtidal zone and others near the sea are the least diversified and structured, while vegetated habitats and surrounding areas are the most diversified and the well-structured.

Relationships between Macrobenthos and Environmental Conditions
Results of the non-parametric multiple regression analysis (DISTLM) between community composition and environmental variables showed significant correlations with nine variables (Adjusted R 2 = 0.41; p < 0.01). These corresponded to the sediment characteristics (mud content, median grain-size, TOM (%), carbon%, CaCO3%, and CO2%), water characteristics (T°, pH, salinity) and habitat type (biomass of Zostera noltei) ( Table 2).  Figure 7 shows the RDA ordination obtained using DISTLM. The pattern indicates that there are at least two trends in the macrofaunal community structure that can be modeled by these environmental drivers. The first clusters, which include stations located near/or in subtidal zones (G6, G7, G8, S4, S10), are driven by salinity, pH, median grain

Relationships between Macrobenthos and Environmental Conditions
Results of the non-parametric multiple regression analysis (DISTLM) between community composition and environmental variables showed significant correlations with nine variables (Adjusted R 2 = 0.41; p < 0.01). These corresponded to the sediment characteristics (mud content, median grain-size, TOM (%), carbon%, CaCO 3 %, and CO 2 %), water characteristics (T • , pH, salinity) and habitat type (biomass of Zostera noltei) ( Table 2).  Figure 7 shows the RDA ordination obtained using DISTLM. The pattern indicates that there are at least two trends in the macrofaunal community structure that can be modeled by these environmental drivers. The first clusters, which include stations located near/or in subtidal zones (G6, G7, G8, S4, S10), are driven by salinity, pH, median grain size and carbon content in the sediment (percentage of: carbon, CO 2 and CaCO 3 ). The second trend highlights the variability between sites in the central and peripheral areas of the lagoon. These variations are related to differences in the percentage of TOM and mud in the sediment, water temperature, salinity and the presence of Z. noltei. The first axis explained 39.8% out of the fitted and 16.3% out of the total variation, while the second accounts for 26.2% of the fitted and 10.7% of the total variation. In total, the first two RDA axes explain 66% of the adjusted change, and this accounts for about 27% of the total change in the multivariate community data. The full RDA axis explains 100% of the adjusted variation and 41.02% of the total variation. size and carbon content in the sediment (percentage of: carbon, CO2 and CaCO3). The second trend highlights the variability between sites in the central and peripheral areas of the lagoon. These variations are related to differences in the percentage of TOM and mud in the sediment, water temperature, salinity and the presence of Z. noltei. The first axis explained 39.8% out of the fitted and 16.3% out of the total variation, while the second accounts for 26.2% of the fitted and 10.7% of the total variation. In total, the first two RDA axes explain 66% of the adjusted change, and this accounts for about 27% of the total change in the multivariate community data. The full RDA axis explains 100% of the adjusted variation and 41.02% of the total variation.

Discussion
Large spatial scale studies are crucial to better manage habitats and resources, particularly for the development of the relatively new ecosystem approach [39][40][41][42]. The relationships between macrobenthos and natural environmental drivers can thus be used to describe habitats, defined as the physical and chemical environment in which a species or community lives, and to provide a baseline for the detection of spatial and temporal changes [43][44][45]. In this context, our study, which covers the entire lagoon, gives an overview on the spatial patterns of the benthic macrofauna of the Moulay Bousselham lagoon in relation to environmental drivers.

Environmental Variables
In this study, not all water parameters varied as gradients between upstream and downstream areas. Only salinity and temperature clearly decreased with gradients upstream in conformity with the findings of [21,46]. Both salinity and temperature gradients resulted from the ocean-continent gradient related to the position of stations across the lagoon and the sampling time.
The salinity of coastal lagoons can vary from freshwater to hypersaline according to local climatic conditions and the degree of hydrological connectivity [47]. However, Figure 7. Two-dimensional redundancy analysis (RDA) ordination representing the model of spatial variation in macrozoobenthos community structure related to the predictor variables selected through the best linear models based on distance (DISTLM).

Discussion
Large spatial scale studies are crucial to better manage habitats and resources, particularly for the development of the relatively new ecosystem approach [39][40][41][42]. The relationships between macrobenthos and natural environmental drivers can thus be used to describe habitats, defined as the physical and chemical environment in which a species or community lives, and to provide a baseline for the detection of spatial and temporal changes [43][44][45]. In this context, our study, which covers the entire lagoon, gives an overview on the spatial patterns of the benthic macrofauna of the Moulay Bousselham lagoon in relation to environmental drivers.

Environmental Variables
In this study, not all water parameters varied as gradients between upstream and downstream areas. Only salinity and temperature clearly decreased with gradients upstream in conformity with the findings of [21,46]. Both salinity and temperature gradients resulted from the ocean-continent gradient related to the position of stations across the lagoon and the sampling time.
The salinity of coastal lagoons can vary from freshwater to hypersaline according to local climatic conditions and the degree of hydrological connectivity [47]. However, within a single lagoon system, there may be three salinity zones whose spatial extent varies depending on seasonal conditions. These are relatively fresh water near the mouths of influent rivers, brackish water in the central part of a lagoon, and marine salinities at the entrance channel (s).
The Moulay Bousselham lagoon is mainly composed of sandy and silty sediments, and no gravel has been detected in each sample measured. The median grain size shows high variations, from fine silts to coarse sands. The main part of the lagoon is composed of poorly sorted sediments according to the Folk and Ward classification, with a mean sorting index of 3.9 estimated from the 68 samples. The sandy dominated stations revealed refer to sediments sampled downstream or in the channels, where higher grain-sized sediments are transported [28], while the lower grain-sized sediments with high mud content identified refer to sediments present in the upstream sections of the lagoon far. They are stations located away from the channels and where morphogenic conditions decrease [48]. Hydrodynamic energy affects sedimentation and resuspension of sediment particles [49,50], as well as organic enrichment of sediments [51,52]. Thus, higher currents and turbulence inhibit the deposition of organic matter and produce the deposition of coarse sediments [49,53], whereas muddy sediments occur in calmer hydrodynamic conditions.
In opposition to the grain size parameters, carbon content presents lower variation, and only low carbon percentages have been recorded by LECO© for the 68 stations studied.
Higher carbon values in stations located nearest channels can be linked to the marine influence providing shell remains in the sands content, whereas the low values (in central and peripheral areas) can be influenced by the continental inputs coming from the watershed and bringing organic matter remains from vegetated areas [29]. The high level of TOM in the center and periphery of the lagoon can be attributed to the presence of fine particles entrapped by the structure of seagrass leaves and the abundance of fragments of dead seagrass encrusted in the sediment [54]. With higher water currents, TOM values are lower near the channels.

Benthic Macrofauna
In the 68 sampled stations, 63 taxa belonging to 50 families were identified. The soft-bottom macrofauna of the Moulay Bousselham lagoon was mainly characterized by the dominance of Mollusca (38.09%), followed by Arthropoda (31.75%), Annelida (28.57%) and Nemertea (1.59%). These results contrast both with previous studies carried out on the lagoon [22,46,55,56], which noted the dominance of mollusks, polychaetes and crustaceans, and with the conclusions of [21], who noted the predominance of crustaceans, followed by polychaetes and mollusks. These results may be related to differences in sampling methods and designs.
On the other hand, the macrobenthic faunal densities observed in this study (4582.8 ind/m 2 ) were higher than those reported by [77] (3106.0 ind/m 2 ) and lower to those observed by [46] (5763.0 ind/m 2 ). The maximum biomass value reported in our study (92.5 g/m 2 ) is higher than that reported by [46]. The highest values were recorded in the vegetated habitats, while the mean value of the biomass (22.2 g/m 2 ) is similar to that reported by [78] (22.0 g/m 2 ) and very close to the results of [46] (20.0 g/m 2 ). Comparisons with these previous studies reveal that the benthic assemblages of the Moulay Bousselham lagoon are relatively stable, indicating a certain durability. Abundance and biomass were clearly lower compared with several other coastal systems: Venice lagoon [72], Prévost lagoon [79], Arcachon bay [80] and the Somme bay [81], but they were higher than those recorded in Boughrara lagoon [82] and Celestun lagoon [67].

Spatial Patterns and Environmental Drivers
The spatial pattern of the benthic communities in the lagoon follows a downstreamupstream gradient, essentially due to environmental factors including sediment characteristics, water parameters and the type of habitat (seagrass beds). The particular combination of those factors generates a macrofaunal structure characterized by 14 assemblages can be clearly seen in the cluster analysis. Assemblages are identified from downstream to upstream and from the center to the peripheral areas. According to cluster analysis, these assemblages showed a clear distinction between the part close to the sea communication (similarity not exceeding 30%) and the parts inside the lagoon (similarity around 60%). Our results also showed the association of stations located in the subtidal with others in the intertidal areas in the identified clusters (G6, G7 and G8). In contrast to the results obtained in the Moulay Bousselham lagoon by [46] and in the Oualidia lagoon by [61], our assemblages do not show a clear dominance of one or two species.
Most of the benthic species inventoried in the Moulay Bousselham lagoon had a wide spatial pattern and were not limited to a single habitat. Such a pattern corresponds better to the concept of a continuum of communities across an environmental gradient [83] than to the concept of discrete communities as distinct assemblages of species defined by [84]. The biological continuum and the absence of ecotonal zones seem to be characteristic of estuaries in particular and semi-enclosed coastal ecosystems [85]. Indeed, this pattern has been found in different estuaries in Morocco [86][87][88], France [89], Portugal [90,91] and Spain [92,93]. The explanation for this finding is probably the high tolerance of the macrozoobenthos species inhabiting such ecosystems. These patterns could also be related to the fact that environmental gradients are not so strong in the Moulay Bousselham lagoon, with the exception of salinity, and that the lagoon lacks large hydrodynamic variations, which commonly have a significant impact on the spatial distribution of benthic communities [94,95]. Analysis of macrobenthic assemblages indicates that the spatial distributions of the 63 species found along the subtidal and the intertidal stations of the Moulay Bousselham lagoon showed a relatively high correlation with environmental drivers and can be best explained by a combination of nine natural abiotic variables. DISTLM highlights sediment characteristics (mud content, median grain-size, TOM%, carbon%, CaCO 3 %, and CO 2 %), water parameters (salinity, T • , pH) and habitat type (biomass of Zostera noltei). There is a gradient from west to east, and the most important stations in terms of specific richness and/or density are those located in the central and peripheral mudflat areas, which are characterized by the presence of a seagrass bed, or located near vegetated areas (Zostera noltei, Ruppia cirrosa, Algae). Past works have shown that the spatial pattern of the benthic communities at Moulay Bousselham lagoon follows an upstream-downstream gradient and demonstrated the primordial role of environmental drivers (sediment grain-size, organic matter, hydrodynamics parameters and the presence of seagrass) on this distribution [21,46].
For macrobenthic invertebrates, such patterns are eventually the result of a complex interaction of a number of processes occurring in both the water column and the sedimentary compartment. Coastal lagoons are complex systems with a high degree of physical and biological variability.
The biodiversity of these ecosystems is commonly thought to be spatially distributed along the vertical and horizontal gradients of salinity, temperature, sediment characteristics (particle size, mud and/or organic matter enrichment) [96][97][98]. This spatial structure results from the environmental tolerances of organisms to stresses within these variable systems (water mass dynamics, physiological stress and biotic interactions) [99,100].
The presence of vegetation creates conditions for the formation of stable and complex habitats, thus promoting the installation of dense and diversified benthic communities [101,102]. The composition of the fauna is also driven by sediment, which is known to be a determinant of macrobenthic composition and plays an important role at different stages of the life cycle (settlement, tube building, burying and feeding) of soft-bottomed benthic organisms [103,104]. Our analyses have also highlighted salinity as a key factor, a parameter that has usually been considered essential to explain gradients in lagoon density, biomass, richness or diversity [105,106] and as one of the main drivers of similarities and differences in lagoon assemblages [107][108][109][110].
Lagoons, however, are characterized by large seasonal, often unpredictable, variation in physical and chemical variables [111,112]. This may act as a driving force regulating the macrozoobenthic assemblages from season to season. At Moulay Bousselham lagoon, previous studies showed that benthic population density and species richness revealed seasonal variation with maxima in the autumn [21]. In the present study, the spatial patterns and associated key environmental drivers were evidenced from sampling performed during autumn where benthic macrofauna are the most diverse. Nevertheless, future studies should consider sampling over different seasons to better trace the physical and biotic factors regulating spatial and seasonal changes in the benthic assemblages of this temperate lagoon.

Conclusions
In conclusion, this study updates the composition of soft-bottom macrofaunal assemblages in the Moulay Bousselham lagoon and provides the first extensive examination of their spatial distribution.
Our results clearly revealed that the hydrographic regime (marine and terrestrial freshwater), the sediment distribution and characteristics, and the type of habitat (vegetated area) are the key factors determining the species composition and patterns of macrozoobenthos assemblages.

Conflicts of Interest:
The authors declare no conflict of interest.