# Latitudinal Variation into the Macrofaunal Assemblages Associated toZostera noltei Seagrass along the Atlantic Coast of Morocco

Authors: Oussama Bououarour; Soilam Boutoumit; Reda El Kamcha; Hocein Bazairi
Journal: International Journal of Zoology
 (2021)

---

## Abstract

Large-scale research on seagrass-associated benthic fauna is very important for future regional marine conservation. In our study, we investigated spatial and latitudinal variation of benthic macroinvertebrate assemblages associated toZostera noltei Hornemann, 1832 beds from five semi-enclosed coastal systems (SECSs) ranging from 23°N to 34°N along the Atlantic coast of Morocco. Overall, 17,320 individuals were reported as belonging to 96 taxa. The ecological community descriptors differ significantly at the level of the site. Specific richness showed an inconsistent significant pattern with latitude. The multivariate analyses of the assemblage’s composition showed 57% of total variation observed in benthic assemblages, while the PERMANOVA analysis confirmed that this variation is significant at the level of the site. According to DistLM results, variations in belowground biomass, and percentage of mud, were the important predictor variables explaining this variation along the large scale of the studied SECS. However, such patterns could be related to other factors such as habitat heterogeneity and regional, biogeographic, and anthropogenic factors. The present study marked the first attempt on broad-scale ecological research of seagrass beds in Morocco and offers baseline data for planning the broad-scale conservation of biodiversity in seagrass beds that remain suffering from multiple human-induced threats such as coastal developments and climate change.


---

## Body

## 1. Introduction

Changes in the composition of community or species assemblages have been and continue to be the subject of intense interest. Anthropogenic pressures and the need for effective conservation planning have further inspired the study of diversity patterns and processes at regional and global scales [1]. Ecologists have been interested in the global pattern of biodiversity for a long time [2], and the comprehension of the distribution of life on the earth is the major goal of ecology and biogeography [3]. The latitudinal diversity gradient (LDG) is one of the most outstanding ecological patterns on our planet. It is generally defined as an increase in species richness from the poles towards the equator and is a striking ecological pattern that has fascinated biologists over centuries [4]. The main drivers behind the LDG are focused on theories that are broadly linked to the current climate, historical effects, and biome area [5, 6].The latitudinal gradient of diversity is well defined in different taxonomic groups and many geographical regions [7, 8], whereas land diversity patterns and their predictors are known for numerous taxa [9], and our understanding of global marine diversity has been more limited, with recent findings revealing some striking contrasts to widely held terrestrial paradigms [10]. Besides, the similarity of species composition is probable to decline with distance because variations in environmental conditions generally upsurge with the distance between regions [11].Seagrass has a wide geographic distribution and presence, except in Antarctica, in many shallow coastal and oceanic waters around the world [12]; therefore, related population and community processes can be compared on a wide spatial scale [13]. They support abundant and generally well-known macrofauna, the density and diversity of which also surpass the existing in nearby bare sediment fields [14, 15].The existence of seagrass meadows in coastal ecosystems favors the establishment and preservation of a high diversity of species in benthic communities [16] which are suitable organisms to test latitudinal patterns as they respond very well to environmental changes and have a high diversity [17]. Seagrass beds and associated invertebrates provide various valuable ecological services including coastal conservation and erosion prevention, carbon sequestration, fisheries maintenance, water purification, and the supply of raw materials and food [18]. However, the quantitative estimation of population and community variables over broad spatial scales is lacking for seagrass-associated communities.The dwarf eelgrassZostera noltei Hornemann, 1832 is the dominant seagrass in the Atlantic semi-enclosed coastal systems (SECSs) of Morocco. It is considered globally as one of the most important perennial seagrass species on intertidal mudflats [19]. It is widely distributed on the Atlantic, from Norway south to the Mauritanian coast, and present in the Mediterranean, the Black Sea, and even in the Caspian and Aral Sea [20]. The existence of Z. noltei with no doubt contributes significantly to the biological, ecological, and environmental values of these SECSs, which are most included in the Ramsar List [21].At the same time, global seagrass beds are declining at a remarkable rate [22]. Zostera noltei is one of the few seagrass species that are adapted to the hard conditions of the intertidal habitat [23]. However, in semi-enclosed coastal areas, this species is most susceptible to impacts resulting from climate change and anthropogenic stresses [24]. Benthic species patterns have rarely been investigated at once on the entirety of the Moroccan coastline (from 20°N to 35°N). Indeed, the existing studies had only a limited geographical scope (Merja Zerga lagoon [25]; Oualidia Lagoon [26]; and Khnifiss Lagoon [27]).In this study, by comparing benthic macrofauna assemblages in five SECSs along the Atlantic coast of Morocco, with consideration of the positions in each SECS (downstream and upstream), we investigated whether there is a latitudinal pattern of benthic macrofauna associated withZ. noltei habitats of Moroccan SECSs. Then, we determined the main environmental factors that could drive to the patterns of distribution and diversity of benthic assemblages. We hypothesized that (1) specific diversity increases towards low latitude along the Moroccan Atlantic coast and (2) that there is a variation in the benthic assemblages among sites and between seagrass positions.


## 2. Materials and Methods

### 2.1. Study Sites

Benthic macrofauna associated toZostera noltei beds were sampled from five SECS distributed along the Atlantic Moroccan coastline: four lagoons: Merja Zerga, Sidi Moussa, Oualidia, and Khnifiss and one bay: Dakhla (Figure 1).Figure 1
Map showing the location of the study sites across the Atlantic coast of Morocco.MZ: Merja Zerga lagoon; SM: Sidi Moussa lagoon; OU: Oualidia lagoon, KH: Khnifiss lagoon, and DA: Dakhla Bay.The Merja Zerga lagoon (34°47′N–6°13′W) is an elliptical-shape lagoon with a depth from 0.50 to 1.50 m. It occupies an area of about 30 km2, and it is largely influenced by tidal rhythms ranging nearly 1.4 m [28]. It is submitted to several pressures such as cattle raising, artisanal fishing, shell-fishing, and tourism [29].The Sidi Moussa lagoon (32°52′N–8°51′W) covers an area of 4.2 km2 with a maximum depth of approximately 5 m, which decreases progressively towards the upstream part of the lagoon. The tidal regime is semidiurnal with a tidal level varying between 2 and 4 m [30]. In the lagoon, various activities such as traditional fishing and aquaculture have a considered impact [31].The Oualidia lagoon (34°47′N–6°13′W) is over 7 km long and 1 km wide, with a mean depth of 2 m for a total surface of 3 km2. It is characterized by a semidiurnal tide and entrances ranging from 0.8 m to 3.6 [32]. The lagoon is subject to various stressors related to fishing activities, aquaculture, and algae exploitation [33].The Khnifiss lagoon (28°02′N–12°13′W) is located at the southern coastal Sahara, and it extends for about 20 km in length, 5 km width for a total surface of 65 km2, and a maximum depth of 8.7 m [34]. The tidal regime is semidiurnal, and it ranges between a minimum of 1.48 m downstream to a maximum of 2.54 m inside the lagoon. In the lagoon, few fishing activities are present [35].The Dakhla Bay (23°35′N–15°50′W) extends over 37 km length and 12 km width for a total surface of 400 km2 and separated from the ocean on its south extremity through a wide 13 km pass [36]. It is a mesotidal system ranging between 0.5 and 2.5 m and with a total depth of no more than 20 m. With the expansion of the harbor, navigation channel dredging, the bay has higher environmental and ecological concerns [37].


### 2.2. Sampling and Analysis

Our sampling was based on a comparison of two seagrass meadows, one downstream and one upstream in each of the five sites between December 2014 and January 2015. In each of the two meadow positions, three stations were randomly sampled (10 replicates per station) in the central dense area of the seagrass bed using a hand PVC corer of 12.5 cm diameter to a depth of 20 cm. The 10 replicates per station totalize a surface of 0.12 m2 per station. Samples were sieved using a mesh of 1 mm and then fixed and conserved with formalin (4%) with Rose Bengal for coloration.Each sample of benthic macrofauna was associated with a sample of sediment collected for the determination of organic carbon and grain-size analyses. The hydrological parameters (water temperature, pH, and salinity) were also measured in situ using a HANNA portable multiparameter. ForZ. noltei meadow characterization, three replicates, using the same hand PVC corer as mentioned above, were sampled and were carefully rinsed on site with seawater to remove remaining inorganic particles and conserved in plastic bags until preparation for analysis.Macrofauna was sorted, identified (to the species level when possible), counted, and preserved in ethanol (70%). The scientific names and the systematic order of species were revised and updated following the WoRMS database (http://marinespecies.org/). Biomass (B) (AFDW, ash free dry weight) was determined following desiccation (48 h at 60°C) and calcination in the oven at 500°C for 3 h.Sediment samples were used for the determination of the different fraction ratios [38]. The grain size was measured using a laser granulometer (Malvern Mastersizer 2000) at LETG, UMR 6554, University of Nantes. Its complete distribution is then treated with the Gradistat© Excel package [39, 40]. To increase the precision of the organic matter estimation, a LECO© carbon analyzer estimates the CO2 and CaCO3 percentages after a 1400°C dioxygen burning and a mineral decarbonizing with a sulfuric acid solution [41, 42].Additionally, three replicates per stationof Zostera noltei beds were randomly sampled for seagrass characterization: shoot density (nb shoots/m2), aboveground biomass (gDW/m2), and belowground biomass (gDW/m2). Samples of Z. noltei were cleaned and separated into leaves and belowground parts (roots and rhizomes) and then were oven-dried at 60°C until constant dry weight (DW).Abundance (N: Ind./m2), number of species (S), Shannon–Wiener diversity index (H′, log2), and Pielou’s evenness index (J′) were calculated for each sample using The DIVERSE routine. A two-way analysis of variance (ANOVA), according to the two-factor (site and position) design, was used to test the differences on S, N, B, H′, and J′ at different scales. Post hoc pairwise multiple comparisons were performed using the Tukey test whenever the interaction between effects showed significant differences (p<0.05). Correlations between these ecological indicators with latitude were tested to determine their pattern along the large gradient.Benthic faunal abundance data were averaged from the three replicates per station. After a fourth root transformation to downweigh the importance of high-abundance species, similarities between sampling stations were calculated using a Bray–Curtis similarity coefficient and then interpreted with the SIMPROF similarity profile test. Environmental variables were log(x+1) transformed and normalized, and a resemblance matrix was created using the Euclidean distance [43]. To visualize differences in overall community structure, Principal Coordinates Ordination (PCO), based on the Bray–Curtis dissimilarity matrix, was performed, as it is considered one of the most suitable visual complements to PERMANOVA output [44], and species that were correlated (Pearson ρ > 0.5) to sample ordination were represented as superimposed vectors in the PCO graph. Differences between sites and positions were tested with a two-way crossed PERMANOVA design, with the site (random factor, with six levels) and position (fixed factor, two levels) being used as factors in the design [45]. The SIMPER routine (cutoff 50%) was used to identify the species most contributing to the similarity of each identified assemblage and the dissimilarity among them [46].Before all statistical analyses, the environmental data were evaluated by draftsman plots to determine collinearity. Tests for collinearity were conducted with no measured collinearity among the environmental parameters (all values < 0.95), and hence, all variables were retained for possible inclusion in the model. We performed a distance-based linear model permutation test (DistLM) to identify which set of environmental variables predicted the multivariate variation in macrofauna assemblages. The adjustedR2 was used as a selection criterion to permit the fitting of the best explanatory environmental variables in the model. Euclidean distance was used as the resemblance measure in all DistLM procedures. Results were visualized using the graphical representation of the distance-based redundancy analysis (dbRDA).All the abovementioned procedures were performed with the PRIMER 6 + PERMANOVA© software (software package from Plymouth Marine Laboratory, UK) [43, 47], while the two-way ANOVA among benthic community structure indexes was carried out in the Statistica software package (StatSoft Inc., 2011, version 10).


## 2.1. Study Sites

Benthic macrofauna associated toZostera noltei beds were sampled from five SECS distributed along the Atlantic Moroccan coastline: four lagoons: Merja Zerga, Sidi Moussa, Oualidia, and Khnifiss and one bay: Dakhla (Figure 1).Figure 1
Map showing the location of the study sites across the Atlantic coast of Morocco.MZ: Merja Zerga lagoon; SM: Sidi Moussa lagoon; OU: Oualidia lagoon, KH: Khnifiss lagoon, and DA: Dakhla Bay.The Merja Zerga lagoon (34°47′N–6°13′W) is an elliptical-shape lagoon with a depth from 0.50 to 1.50 m. It occupies an area of about 30 km2, and it is largely influenced by tidal rhythms ranging nearly 1.4 m [28]. It is submitted to several pressures such as cattle raising, artisanal fishing, shell-fishing, and tourism [29].The Sidi Moussa lagoon (32°52′N–8°51′W) covers an area of 4.2 km2 with a maximum depth of approximately 5 m, which decreases progressively towards the upstream part of the lagoon. The tidal regime is semidiurnal with a tidal level varying between 2 and 4 m [30]. In the lagoon, various activities such as traditional fishing and aquaculture have a considered impact [31].The Oualidia lagoon (34°47′N–6°13′W) is over 7 km long and 1 km wide, with a mean depth of 2 m for a total surface of 3 km2. It is characterized by a semidiurnal tide and entrances ranging from 0.8 m to 3.6 [32]. The lagoon is subject to various stressors related to fishing activities, aquaculture, and algae exploitation [33].The Khnifiss lagoon (28°02′N–12°13′W) is located at the southern coastal Sahara, and it extends for about 20 km in length, 5 km width for a total surface of 65 km2, and a maximum depth of 8.7 m [34]. The tidal regime is semidiurnal, and it ranges between a minimum of 1.48 m downstream to a maximum of 2.54 m inside the lagoon. In the lagoon, few fishing activities are present [35].The Dakhla Bay (23°35′N–15°50′W) extends over 37 km length and 12 km width for a total surface of 400 km2 and separated from the ocean on its south extremity through a wide 13 km pass [36]. It is a mesotidal system ranging between 0.5 and 2.5 m and with a total depth of no more than 20 m. With the expansion of the harbor, navigation channel dredging, the bay has higher environmental and ecological concerns [37].


## 2.2. Sampling and Analysis

Our sampling was based on a comparison of two seagrass meadows, one downstream and one upstream in each of the five sites between December 2014 and January 2015. In each of the two meadow positions, three stations were randomly sampled (10 replicates per station) in the central dense area of the seagrass bed using a hand PVC corer of 12.5 cm diameter to a depth of 20 cm. The 10 replicates per station totalize a surface of 0.12 m2 per station. Samples were sieved using a mesh of 1 mm and then fixed and conserved with formalin (4%) with Rose Bengal for coloration.Each sample of benthic macrofauna was associated with a sample of sediment collected for the determination of organic carbon and grain-size analyses. The hydrological parameters (water temperature, pH, and salinity) were also measured in situ using a HANNA portable multiparameter. ForZ. noltei meadow characterization, three replicates, using the same hand PVC corer as mentioned above, were sampled and were carefully rinsed on site with seawater to remove remaining inorganic particles and conserved in plastic bags until preparation for analysis.Macrofauna was sorted, identified (to the species level when possible), counted, and preserved in ethanol (70%). The scientific names and the systematic order of species were revised and updated following the WoRMS database (http://marinespecies.org/). Biomass (B) (AFDW, ash free dry weight) was determined following desiccation (48 h at 60°C) and calcination in the oven at 500°C for 3 h.Sediment samples were used for the determination of the different fraction ratios [38]. The grain size was measured using a laser granulometer (Malvern Mastersizer 2000) at LETG, UMR 6554, University of Nantes. Its complete distribution is then treated with the Gradistat© Excel package [39, 40]. To increase the precision of the organic matter estimation, a LECO© carbon analyzer estimates the CO2 and CaCO3 percentages after a 1400°C dioxygen burning and a mineral decarbonizing with a sulfuric acid solution [41, 42].Additionally, three replicates per stationof Zostera noltei beds were randomly sampled for seagrass characterization: shoot density (nb shoots/m2), aboveground biomass (gDW/m2), and belowground biomass (gDW/m2). Samples of Z. noltei were cleaned and separated into leaves and belowground parts (roots and rhizomes) and then were oven-dried at 60°C until constant dry weight (DW).Abundance (N: Ind./m2), number of species (S), Shannon–Wiener diversity index (H′, log2), and Pielou’s evenness index (J′) were calculated for each sample using The DIVERSE routine. A two-way analysis of variance (ANOVA), according to the two-factor (site and position) design, was used to test the differences on S, N, B, H′, and J′ at different scales. Post hoc pairwise multiple comparisons were performed using the Tukey test whenever the interaction between effects showed significant differences (p<0.05). Correlations between these ecological indicators with latitude were tested to determine their pattern along the large gradient.Benthic faunal abundance data were averaged from the three replicates per station. After a fourth root transformation to downweigh the importance of high-abundance species, similarities between sampling stations were calculated using a Bray–Curtis similarity coefficient and then interpreted with the SIMPROF similarity profile test. Environmental variables were log(x+1) transformed and normalized, and a resemblance matrix was created using the Euclidean distance [43]. To visualize differences in overall community structure, Principal Coordinates Ordination (PCO), based on the Bray–Curtis dissimilarity matrix, was performed, as it is considered one of the most suitable visual complements to PERMANOVA output [44], and species that were correlated (Pearson ρ > 0.5) to sample ordination were represented as superimposed vectors in the PCO graph. Differences between sites and positions were tested with a two-way crossed PERMANOVA design, with the site (random factor, with six levels) and position (fixed factor, two levels) being used as factors in the design [45]. The SIMPER routine (cutoff 50%) was used to identify the species most contributing to the similarity of each identified assemblage and the dissimilarity among them [46].Before all statistical analyses, the environmental data were evaluated by draftsman plots to determine collinearity. Tests for collinearity were conducted with no measured collinearity among the environmental parameters (all values < 0.95), and hence, all variables were retained for possible inclusion in the model. We performed a distance-based linear model permutation test (DistLM) to identify which set of environmental variables predicted the multivariate variation in macrofauna assemblages. The adjustedR2 was used as a selection criterion to permit the fitting of the best explanatory environmental variables in the model. Euclidean distance was used as the resemblance measure in all DistLM procedures. Results were visualized using the graphical representation of the distance-based redundancy analysis (dbRDA).All the abovementioned procedures were performed with the PRIMER 6 + PERMANOVA© software (software package from Plymouth Marine Laboratory, UK) [43, 47], while the two-way ANOVA among benthic community structure indexes was carried out in the Statistica software package (StatSoft Inc., 2011, version 10).


## 3. Results

### 3.1. Environmental Variables

Values of environmental parameters are shown in Table1. The mean temperature is ranging from 14.06 to 21.38°C with a clear variation between upstream and downstream samples. pH values were different between sites in both positions with a variation among position areas. Salinity values present a significant variation beyond sites ranging from 16.25 to 39.75. The sediment composition revealed that the mud content is varying from 8.00 to 88.43% in the sampling sites. Concentration of CO2 and CaCO3 showed a significant difference among sites (p<0.05) with no difference between downstream and upstream samples.Table 1
Average values of environmental variables recorded among the study sites from two sampling position habitats.
SiteMerja Zerga lagoonSidi Moussa lagoonOualidia lagoonKhnifiss lagoonDakhla BayPositionDownstreamUpstreamDownstreamUpstreamDownstreamUpstreamDownstreamUpstreamDownstreamUpstreamLatitude (N)34.8734.8532.9732.9932.7432.7628.0728.0823.8123.66Water parametersTemperature (°C)14.0616.2215.5817.3815.4816.3121.3817.0920.7416.93pH8.919.438.778.038.228.608.797.817.848.25Salinity16.2518.8629.6624.3434.2229.4836.2639.7539.4427.02Sediment parametersMud fraction (%)79.8170.7268.0188.4315.0366.4665.3081.808.0019.70Sand fraction (%)20.1929.2831.9811.5784.9633.5334.7018.2092.0080.30CO2 (%)35.53 ± 0.3134.63 ± 2.4875.4 ± 4.1346.53 ± 6.5599.43 ± 0.3164.46 ± 8.6250.83 ± 21.7789.43 ± 13.2246.03 ± 7.0417.43 ± 1.17CaCO3 (%)80.86 ± 0.6878.76 ± 5.5798.13 ± 0.2292.4 ± 5.9399.43 ± 0.3197.63 ± 0.3785.1 ± 8.9398.83 ± 0.6894.16 ± 2.1739.73 ± 2.71Seagrass metricsShoot density (shoot. m−2)2194.33 ± 184.83388.6 ± 425.72944 ± 907.35027.3 ± 1592.44929.6 ± 1380.22944 ± 1240.63462.7 ± 253.14620.2 ± 283.95425.8 ± 2938.26249.8 ± 444.2Aboveground biomass (g DW m−2)67.8 ± 22.775.7 ± 6.394.3 ± 17.6146.1 ± 39.294.5 ± 43.787.1 ± 21.981.6 ± 15.456.4 ± 7.187.5 ± 36.777.1 ± 7.3Belowground biomass (g DW m−2)46.8 ± 11.143.4 ± 16.738.5 ± 16.134.3 ± 11.222.8 ± 3.331.3 ± 2.296.9 ± 8.586.9 ± 6.76233.8 ± 53.8223.9 ± 29.8For seagrass measurements, shoot density ranged from 2194 to 6250 shoots. m−2, while aboveground biomass was varying between 56.4 and 146.1 g DW m−2 and the belowground biomass was fluctuating from 22.8 to 223.9 g DW m−2. All the seagrass parameters did not present a significant difference at the level of position (downstream/upstream) (p>0.05), while at the level of sites, only the belowground biomass was significantly different (p<0.05). The Pearson correlation between latitude and all the abovementioned environmental variables revealed a significant relationship between temperature (r = −0.672), shoot density (r = −0.716), and belowground biomass (r = −0.931) (Table 2).Table 2
Pearson’s correlation coefficients (r) for the relationships between latitude and each of the environmental variables used for this study.
Correlation with latitude (N)Temperature (°C)−0.672pH0.598Salinity−0.631Mud (%)0.601CO2 (%)0.222CaCO3 (%)0.404Shoot density (shoot. m−2)−0.716Aboveground biomass (g DW m−2)0.218Belowground biomass (g DW m−2)−0.931Bold values mean that significative correlations based onp < 0.05 were detected for the tested parameter.


### 3.2. Diversity and Species Composition

This study identified a total of 17,320 individuals belonging to 96 benthic macroinvertebrate taxa. At the level of the phylum, theZostera noltei beds in the sampled SECS were highly dominated by crustaceans with 38 species, representing more than 38% of the sampled macrofauna. The mollusks were composed of 24 species representing 33% of the total of individuals followed by the polychaetes with nearly 25% of the abundance belonging to 27 species. The other phyla (platyhelminths, nemerteans, echinoderms, and cnidarians) were far less abundant with 3.18% of the mean global abundance with only seven species.Taxonomic richness reached the maximum values in the Dakhla downstream sample (DA-D) (31 ± 3), and the lowest mean abundance was recorded in the Merja Zerga downstream site (MZ-D) (11 ± 1). Density ranged from 2361 ± 264 individuals per m2 in Oualidia upstream samples (site OU-U) to 9464 ± 902 individuals per m2 in the Merja Zerga downstream site (MZ-D). Biomass was greater in Merja Zerga upstream (MZ-U) with 244.5 ± 48.7 g AFDW m−2, while the lowest biomass values were noted in Khnifiss downstream samples (KH-D) (1.31 ± 0.1 g AFDW m−2) (Table 3).Table 3
Mean values of ecological indices for samples occuping two positions (downstream and upstream) in the five study sites.
SitesMerja Zerga lagoonSidi Moussa lagoonOualidia lagoonKhnifiss lagoonDakhla BayPositionDownstreamUpstreamDownstreamUpstreamDownstreamUpstreamDownstreamUpstreamDownstreamUpstreamS11 ± 113 ± 322 ± 317 ± 220 ± 114 ± 215 ± 320 ± 125 ± 231 ± 3N9464 ± 9024244 ± 4982427.77 ± 8092478 ± 7595039 ± 4312361 ± 2463389 ± 431.484625 ± 8837264 ± 2986819 ± 418B210.2 ± 158.55244.54 ± 48.7216.06 ± 11.2490.83 ± 25.2160.81 ± 35.2817.65 ± 6.771.31 ± 0.173.49 ± 0.5245.46 ± 1.0211.98 ± 3.01H′1.179 ± 0.281.99 ± 0.693.33 ± 0.222.81 ± 0.293.13 ± 0.041.88 ± 0.112.55 ± 0.193.18 ± 0.152.76 ± 0.223.87 ± 0.12J′0.34 ± 0.070.52 ± 0.140.74 ± 0.050.69 ± 0.070.72 ± 0.020.49 ± 0.050.64 ± 0.020.73 ± 0.030.59 ± 0.030.78 ± 0.01S- species richness per 0.36 m², N- abundance per m2, B- biomass per g AFDW m−2, H'– Shannon–Wiener diversity index (log2), J′- Pielou`s evenness index.Taxonomic richness (S), abundance (N), biomass (B), Shannon’s diversity index (H′), and Pielou’s evenness (J′) were significantly different at the level of sites (Table 4). Post hoc analyses showed that all these descriptors’ values are increasing significantly from Merja Zerga lagoon to Dakhla Bay. However, at the level of position, the only descriptor that did not show a significant difference was the biomass. Interactions between the two effects showed significant differences for N, H′, and J′ (Table 4). The mean species richness of benthic macrofauna from all the stations investigated in the present study decreased with latitude (R = −0.76; p=0.01) (Figure 2). The other ecological descriptors did not show a significant correlation with latitude.Table 4
Results of two-way ANOVA testing for macrobenthic assemblages’ differences in the total number of individuals (N), total number of taxa (S), biomass (B), Shannon’s diversity index (H′), and Pielou’s evenness index (J′) at the scales of site and position (downstream-upstream) with post hoc pairwise multiple comparisons (Tukey test) among the studied sites.
SNBH′J′DfFpFpFpFpFpSite422.330.00136.690.0019.450.00116.970.00110.930.001Position113.480.00133.010.0010.550.46513.410.0015.390.031Site × Position42.420.08111.720.0010.570.6826.220.0025.010.005Error20Total29Tukey’s HSD for site × positionIn terms of “Position” effectDMZ ≠ SM ≠ OU = KH = DAMZ ≠ SM = OU = KH = DAMZ ≠ SM = OU = KH = DAUMZ = SM = OU = KH ≠DAMZ = SM = OU = KH = DAMZ = SM = OU = KH = DABold values mean that significative differences based onp < 0.05 were detected for the tested effect.Figure 2
Pearson’s correlation analyses of mean species richness ofZ. noltei beds present in the sampled sites with latitude. Samples were collected from downstream and upstream stations of five semi-enclosed coastal ecosystems along the Atlantic coast of Morocco. Grey shaded area: 95% confidence intervals of best-fit lines.


### 3.3. Latitudinal Patterns of Benthic Assemblages

PERMANOVA analysis showed significant and independent differences in the structure of macrofauna for site (p perm < 0.05), while there was no significant difference in terms of position (p perm > 0.05). The interaction between effects was significant (p perm = 0.0001) (Table 5).Table 5
Results of the multivariate permutational analysis (PERMANOVA) for benthic assemblages considering “site” (5 levels, random) and “position” (2 levels, fixed) factors. Analyses are based on Bray–Curtis similarity matrixes from fourth root transformed data; Df = degrees of freedom; MS = mean square; andp (perm) = level of significance.
SourcedfMSPseudo-Fp (perm)Site41080918.7020.001Position15246.41.16110.367Site x position44518.37.81780.001Residual20577.94Total29Pairwise comparisons for the term site x positionFor pairs of levels of “site” within the level of “position”DMZ ≠ SM = OU ≠ KH ≠ DAUMZ ≠ SM ≠ OU ≠ KH ≠ DABold values mean that significative differences based on p < 0.05 were detected for the tested variable.SIMPER analysis showed that the species contributing the most to the average similarity within group I (60.51%) were the two mollusksPeringia ulvae Pennant, 1777 and Scrobicularia plana da Costa, 1778, the polychaete Heteromastus filiformis Claparède, 1864, and the isopod Idotea chelipes Pallas, 1766; group II (56.36% of average similarity) is characterized by the presence of the gastropod P. ulvae, three bivalves (Cerastoderma edule Linnaeus, 1758, S. plana, and Abra tenuis Montagu, 1803), the crustaceans I. chelipes and Microdeutopus chelifer Spence Bate, 1862, the polychaete H. filiformis, and the cnidarian Actinia equina Linnaeus, 1758. Group III (63.52% of average similarity) is dominated mainly by crustacean species, two isopods (Cyathura carinata Krøyer, 1847 and I. chelipes), and two amphipods (Lysianassa ceratina Walker, 1889 and Melita palmata Montagu, 1804) with the contribution of the polychaetes (Cirriformia tentaculata Montagu, 1808, Aphelochaeta sp. Blake, 1991, and Diopatra sp. Audouin and Milne Edwards, 1833). Finally, group IV was defined by four crustaceans (C. carinata, I. chelipes, M. chelifer, and Urothoe elegans Spence Bate, 1857) and five polychaetes (Euclymene palermitana Grube, 1840, H. diversicolor, Scoloplos armiger Müller, 1776, Nicomache (Loxochona) trispinata Arwidsson, 1906, and H. filiformis) (Table 6).Table 6
Results of SIMPER analysis showing the average similarity within the different benthic assemblages. Assemblages were identified by the hierarchical ascendant classification analysis and the characteristic species of each benthic assemblage.
GroupsG IG IIG IIIG IVAverage similarity (%)60.5156.3663.5253.11Species contribution, %Abra tenuis6.89Actinia equina6.60Aphelochaeta sp.6.53Cerastoderma edule7.09Cirriformia tentaculata8.35Cyathura carinata10.485.64Diopatra sp.6.35Euclymene palermitana7.16Hediste diversicolor5.79Heteromastus filiformis13.082.755.10Idotea chelipes10.707.647.235.37Lysianassa ceratina7.97Melita palmata6.41Microdeutopus chelifer1.915.61Nicomache (Loxochona) trispinata5.37Peringia ulvae17.608.83Scoloplos armiger5.48Scrobicularia plana13.776.96Urothoe elegans5.44Bold values refer to the highest contributive species on the similarity in each benthic assemblage.The dissimilarity between group I and the other groups (G II, 56.86%; G III, 79.33%, and G IV 76.03%) was mainly due to the dominance of the gastropodP. ulvae and the amphipod M. chelifer. The main responsible taxa of dissimilarity of group II with other groups were A. tenuis, C. edule, C. carinata, M. palmata, and Ampithoe ferox Chevreux, 1901. The species Leucothoe richiardii Lesson, 1865, Aphelochaeta sp., and Terebella lapidaria Linnaeus, 1767 are the main contributors in the dissimilarity of group III with the other groups. Likewise, the dissimilarity-contributing taxa between group IV and the other assemblages were the crustaceans Ampelisca sp Krøyer, 1842, Apseudes sp. Leach, 1814, Caprella acanthifera Leach, 1814, and U. elegans, with N. trispinata and S. armiger (Table 7).Table 7
Results of SIMPER analysis showing the average dissimilarity between the different benthic assemblages. Assemblages were identified by the hierarchical ascendant classification analysis and the species contributing in the dissimilarity of each benthic assemblage (cutoff 50%).
GroupsG I and G IIG I and G IIIG I and G IVG II and G IIIG II and G IVG III and G IVAverage dissimilarity (%)56.8673.2276.0367.0472.7166.01Species contribution, %Abra tenuis3.603.872.78Actinia equina4.233.092.362.60Ampelisca sp.3.022.863.12Amphitoe ramondi2.782.62Ampithoe ferox3.762.782.03Aphelochaeta sp.3.813.662.87Apseudes latreillii1.51Apseudes sp.2.542.422.64Bittium reticulatum2.15Caprella acanthifera3.493.273.57Caprella takeuchii1.651.80Cerastoderma edule3.122.491.73Chone duneri1.911.571.67Cirriformia tentaculata4.384.211.54Clibanarius erythropus2.49Cyathura carinata3.884.782.23Diopatra sp.4.133.233.973.06Dynamene sp.1.53Euclymene palermitana3.133.592.321.713.71Hediste diversicolorHeteromastus filiformis3.705.192.262.40Lysianassa ceratina3.054.651.623.51Lekanesphaera rugicauda3.561.98Leucothoe richiardii2.922.792.17Loripes orbiculatus1.85Lagis koreni2.69Melita palmata2.582.231.641.68Microdeutopus chelifer4.513.623.381.56Nicomache (Loxochona) trispinata2.722.572.80Ophiura sp.2.951.81Peringia ulvae4.468.396.104.983.55Sabella spp.2.431.91Scoloplos armiger2.362.353.06Scrobicularia plana2.803.812.561.84Tanais dulongii2.722.05Terebella lapidaria5.184.952.87Tritia pfeifferi1.68Urothoe elegans3.062.893.16The first two axes of the PCO analyses explained 57% of total variation, and the PCO plot indicated a distinct pattern of the benthic communities’ structure between the “northern” and the “southern” site’s samples (Figure3). They were clearly separated along the PCO1 axis (37.4%) which is negatively correlated with C. tentaculata (r = −0.83), Chone duneri Malmgren, 1867 (r = −0.72), Ampithoe ramondi Audouin, 1826 (r = −0.91), and Ampelisca sp. (r = −0.64) and correlated positively with A. tenuis (0.68), while the PCO2 axis (19.6% of total of variation) separating between groups in both regions is negatively correlated positively with Bittium reticulatum da Costa, 1778 (r = 0.61), Capitella capitata Fabricius, 1780 (r = 0.66), Ampelisca sp. (r = 0.67), Apseudes sp. (r = 0.63), and negatively with Aphelochaeta sp. (r = −0.72).Figure 3
Principal coordinate analysis (PCO) of benthic fauna mean abundance data of sampled sites. Attributes that were most associated to the PCO axes (Person correlations > 0.3) are represented by vectors that overlap the PCO.


### 3.4. Relationship between Environmental and Biological Data

The sequential DistLM analysis showed that the belowground biomass of the seagrass, the salinity, and the percentage of mud fraction in the sediment had a significant correlation on the latitudinal distribution of benthic assemblages (p<0.05), explaining the greatest proportion (63%). However, the best solution provided through the DistLM analysis was found when using five variables (salinity, temperature, percentage of mud, aboveground, and belowground biomasses) as environmental predictors of benthic macrofauna composition, explaining 78% of the total variability between samples (Table 8). The first two dbRDA axes captured 69.8% of the variability in the fitted model and 54.6% of the total variation in the data cloud (Figure 4). The first dbRDA1 axis (36% of the total variation) is correlated with belowground biomass (r = −0.71) and % salinity (r = −0.46). The dbRDA2 axis represents 18.6% of the total variation and correlated strongly with the percentage of mud (r = −0.61) and salinity (r = −0.6) (Figure 4).Table 8
Results of the sequential test of the multivariate regression analysis (DistLM). Environmental variables were analyzed individually (marginal test) and sequentially using a forward selection procedure (adjustedR2 = an information criterion).
VariablePseudo-Fp valueProportionCumulativeproportionBelowground biomass3.83360.00040.32390.3237Salinity2.64180.00830.18520.5092Percentage of mud2.0680.04740.12580.6350Aboveground biomass1.26970.28820.07390.7089Temperature1.36130.26050.07390.7828Figure 4
Two-dimensional distance-based redundancy analysis (dbRDA) ordination representing the model of spatial variation in macrofaunal community structure related to the predictor variables selected through the best linear models based on distance (DistLM).


## 3.1. Environmental Variables

Values of environmental parameters are shown in Table1. The mean temperature is ranging from 14.06 to 21.38°C with a clear variation between upstream and downstream samples. pH values were different between sites in both positions with a variation among position areas. Salinity values present a significant variation beyond sites ranging from 16.25 to 39.75. The sediment composition revealed that the mud content is varying from 8.00 to 88.43% in the sampling sites. Concentration of CO2 and CaCO3 showed a significant difference among sites (p<0.05) with no difference between downstream and upstream samples.Table 1
Average values of environmental variables recorded among the study sites from two sampling position habitats.
SiteMerja Zerga lagoonSidi Moussa lagoonOualidia lagoonKhnifiss lagoonDakhla BayPositionDownstreamUpstreamDownstreamUpstreamDownstreamUpstreamDownstreamUpstreamDownstreamUpstreamLatitude (N)34.8734.8532.9732.9932.7432.7628.0728.0823.8123.66Water parametersTemperature (°C)14.0616.2215.5817.3815.4816.3121.3817.0920.7416.93pH8.919.438.778.038.228.608.797.817.848.25Salinity16.2518.8629.6624.3434.2229.4836.2639.7539.4427.02Sediment parametersMud fraction (%)79.8170.7268.0188.4315.0366.4665.3081.808.0019.70Sand fraction (%)20.1929.2831.9811.5784.9633.5334.7018.2092.0080.30CO2 (%)35.53 ± 0.3134.63 ± 2.4875.4 ± 4.1346.53 ± 6.5599.43 ± 0.3164.46 ± 8.6250.83 ± 21.7789.43 ± 13.2246.03 ± 7.0417.43 ± 1.17CaCO3 (%)80.86 ± 0.6878.76 ± 5.5798.13 ± 0.2292.4 ± 5.9399.43 ± 0.3197.63 ± 0.3785.1 ± 8.9398.83 ± 0.6894.16 ± 2.1739.73 ± 2.71Seagrass metricsShoot density (shoot. m−2)2194.33 ± 184.83388.6 ± 425.72944 ± 907.35027.3 ± 1592.44929.6 ± 1380.22944 ± 1240.63462.7 ± 253.14620.2 ± 283.95425.8 ± 2938.26249.8 ± 444.2Aboveground biomass (g DW m−2)67.8 ± 22.775.7 ± 6.394.3 ± 17.6146.1 ± 39.294.5 ± 43.787.1 ± 21.981.6 ± 15.456.4 ± 7.187.5 ± 36.777.1 ± 7.3Belowground biomass (g DW m−2)46.8 ± 11.143.4 ± 16.738.5 ± 16.134.3 ± 11.222.8 ± 3.331.3 ± 2.296.9 ± 8.586.9 ± 6.76233.8 ± 53.8223.9 ± 29.8For seagrass measurements, shoot density ranged from 2194 to 6250 shoots. m−2, while aboveground biomass was varying between 56.4 and 146.1 g DW m−2 and the belowground biomass was fluctuating from 22.8 to 223.9 g DW m−2. All the seagrass parameters did not present a significant difference at the level of position (downstream/upstream) (p>0.05), while at the level of sites, only the belowground biomass was significantly different (p<0.05). The Pearson correlation between latitude and all the abovementioned environmental variables revealed a significant relationship between temperature (r = −0.672), shoot density (r = −0.716), and belowground biomass (r = −0.931) (Table 2).Table 2
Pearson’s correlation coefficients (r) for the relationships between latitude and each of the environmental variables used for this study.
Correlation with latitude (N)Temperature (°C)−0.672pH0.598Salinity−0.631Mud (%)0.601CO2 (%)0.222CaCO3 (%)0.404Shoot density (shoot. m−2)−0.716Aboveground biomass (g DW m−2)0.218Belowground biomass (g DW m−2)−0.931Bold values mean that significative correlations based onp < 0.05 were detected for the tested parameter.


## 3.2. Diversity and Species Composition

This study identified a total of 17,320 individuals belonging to 96 benthic macroinvertebrate taxa. At the level of the phylum, theZostera noltei beds in the sampled SECS were highly dominated by crustaceans with 38 species, representing more than 38% of the sampled macrofauna. The mollusks were composed of 24 species representing 33% of the total of individuals followed by the polychaetes with nearly 25% of the abundance belonging to 27 species. The other phyla (platyhelminths, nemerteans, echinoderms, and cnidarians) were far less abundant with 3.18% of the mean global abundance with only seven species.Taxonomic richness reached the maximum values in the Dakhla downstream sample (DA-D) (31 ± 3), and the lowest mean abundance was recorded in the Merja Zerga downstream site (MZ-D) (11 ± 1). Density ranged from 2361 ± 264 individuals per m2 in Oualidia upstream samples (site OU-U) to 9464 ± 902 individuals per m2 in the Merja Zerga downstream site (MZ-D). Biomass was greater in Merja Zerga upstream (MZ-U) with 244.5 ± 48.7 g AFDW m−2, while the lowest biomass values were noted in Khnifiss downstream samples (KH-D) (1.31 ± 0.1 g AFDW m−2) (Table 3).Table 3
Mean values of ecological indices for samples occuping two positions (downstream and upstream) in the five study sites.
SitesMerja Zerga lagoonSidi Moussa lagoonOualidia lagoonKhnifiss lagoonDakhla BayPositionDownstreamUpstreamDownstreamUpstreamDownstreamUpstreamDownstreamUpstreamDownstreamUpstreamS11 ± 113 ± 322 ± 317 ± 220 ± 114 ± 215 ± 320 ± 125 ± 231 ± 3N9464 ± 9024244 ± 4982427.77 ± 8092478 ± 7595039 ± 4312361 ± 2463389 ± 431.484625 ± 8837264 ± 2986819 ± 418B210.2 ± 158.55244.54 ± 48.7216.06 ± 11.2490.83 ± 25.2160.81 ± 35.2817.65 ± 6.771.31 ± 0.173.49 ± 0.5245.46 ± 1.0211.98 ± 3.01H′1.179 ± 0.281.99 ± 0.693.33 ± 0.222.81 ± 0.293.13 ± 0.041.88 ± 0.112.55 ± 0.193.18 ± 0.152.76 ± 0.223.87 ± 0.12J′0.34 ± 0.070.52 ± 0.140.74 ± 0.050.69 ± 0.070.72 ± 0.020.49 ± 0.050.64 ± 0.020.73 ± 0.030.59 ± 0.030.78 ± 0.01S- species richness per 0.36 m², N- abundance per m2, B- biomass per g AFDW m−2, H'– Shannon–Wiener diversity index (log2), J′- Pielou`s evenness index.Taxonomic richness (S), abundance (N), biomass (B), Shannon’s diversity index (H′), and Pielou’s evenness (J′) were significantly different at the level of sites (Table 4). Post hoc analyses showed that all these descriptors’ values are increasing significantly from Merja Zerga lagoon to Dakhla Bay. However, at the level of position, the only descriptor that did not show a significant difference was the biomass. Interactions between the two effects showed significant differences for N, H′, and J′ (Table 4). The mean species richness of benthic macrofauna from all the stations investigated in the present study decreased with latitude (R = −0.76; p=0.01) (Figure 2). The other ecological descriptors did not show a significant correlation with latitude.Table 4
Results of two-way ANOVA testing for macrobenthic assemblages’ differences in the total number of individuals (N), total number of taxa (S), biomass (B), Shannon’s diversity index (H′), and Pielou’s evenness index (J′) at the scales of site and position (downstream-upstream) with post hoc pairwise multiple comparisons (Tukey test) among the studied sites.
SNBH′J′DfFpFpFpFpFpSite422.330.00136.690.0019.450.00116.970.00110.930.001Position113.480.00133.010.0010.550.46513.410.0015.390.031Site × Position42.420.08111.720.0010.570.6826.220.0025.010.005Error20Total29Tukey’s HSD for site × positionIn terms of “Position” effectDMZ ≠ SM ≠ OU = KH = DAMZ ≠ SM = OU = KH = DAMZ ≠ SM = OU = KH = DAUMZ = SM = OU = KH ≠DAMZ = SM = OU = KH = DAMZ = SM = OU = KH = DABold values mean that significative differences based onp < 0.05 were detected for the tested effect.Figure 2
Pearson’s correlation analyses of mean species richness ofZ. noltei beds present in the sampled sites with latitude. Samples were collected from downstream and upstream stations of five semi-enclosed coastal ecosystems along the Atlantic coast of Morocco. Grey shaded area: 95% confidence intervals of best-fit lines.


## 3.3. Latitudinal Patterns of Benthic Assemblages

PERMANOVA analysis showed significant and independent differences in the structure of macrofauna for site (p perm < 0.05), while there was no significant difference in terms of position (p perm > 0.05). The interaction between effects was significant (p perm = 0.0001) (Table 5).Table 5
Results of the multivariate permutational analysis (PERMANOVA) for benthic assemblages considering “site” (5 levels, random) and “position” (2 levels, fixed) factors. Analyses are based on Bray–Curtis similarity matrixes from fourth root transformed data; Df = degrees of freedom; MS = mean square; andp (perm) = level of significance.
SourcedfMSPseudo-Fp (perm)Site41080918.7020.001Position15246.41.16110.367Site x position44518.37.81780.001Residual20577.94Total29Pairwise comparisons for the term site x positionFor pairs of levels of “site” within the level of “position”DMZ ≠ SM = OU ≠ KH ≠ DAUMZ ≠ SM ≠ OU ≠ KH ≠ DABold values mean that significative differences based on p < 0.05 were detected for the tested variable.SIMPER analysis showed that the species contributing the most to the average similarity within group I (60.51%) were the two mollusksPeringia ulvae Pennant, 1777 and Scrobicularia plana da Costa, 1778, the polychaete Heteromastus filiformis Claparède, 1864, and the isopod Idotea chelipes Pallas, 1766; group II (56.36% of average similarity) is characterized by the presence of the gastropod P. ulvae, three bivalves (Cerastoderma edule Linnaeus, 1758, S. plana, and Abra tenuis Montagu, 1803), the crustaceans I. chelipes and Microdeutopus chelifer Spence Bate, 1862, the polychaete H. filiformis, and the cnidarian Actinia equina Linnaeus, 1758. Group III (63.52% of average similarity) is dominated mainly by crustacean species, two isopods (Cyathura carinata Krøyer, 1847 and I. chelipes), and two amphipods (Lysianassa ceratina Walker, 1889 and Melita palmata Montagu, 1804) with the contribution of the polychaetes (Cirriformia tentaculata Montagu, 1808, Aphelochaeta sp. Blake, 1991, and Diopatra sp. Audouin and Milne Edwards, 1833). Finally, group IV was defined by four crustaceans (C. carinata, I. chelipes, M. chelifer, and Urothoe elegans Spence Bate, 1857) and five polychaetes (Euclymene palermitana Grube, 1840, H. diversicolor, Scoloplos armiger Müller, 1776, Nicomache (Loxochona) trispinata Arwidsson, 1906, and H. filiformis) (Table 6).Table 6
Results of SIMPER analysis showing the average similarity within the different benthic assemblages. Assemblages were identified by the hierarchical ascendant classification analysis and the characteristic species of each benthic assemblage.
GroupsG IG IIG IIIG IVAverage similarity (%)60.5156.3663.5253.11Species contribution, %Abra tenuis6.89Actinia equina6.60Aphelochaeta sp.6.53Cerastoderma edule7.09Cirriformia tentaculata8.35Cyathura carinata10.485.64Diopatra sp.6.35Euclymene palermitana7.16Hediste diversicolor5.79Heteromastus filiformis13.082.755.10Idotea chelipes10.707.647.235.37Lysianassa ceratina7.97Melita palmata6.41Microdeutopus chelifer1.915.61Nicomache (Loxochona) trispinata5.37Peringia ulvae17.608.83Scoloplos armiger5.48Scrobicularia plana13.776.96Urothoe elegans5.44Bold values refer to the highest contributive species on the similarity in each benthic assemblage.The dissimilarity between group I and the other groups (G II, 56.86%; G III, 79.33%, and G IV 76.03%) was mainly due to the dominance of the gastropodP. ulvae and the amphipod M. chelifer. The main responsible taxa of dissimilarity of group II with other groups were A. tenuis, C. edule, C. carinata, M. palmata, and Ampithoe ferox Chevreux, 1901. The species Leucothoe richiardii Lesson, 1865, Aphelochaeta sp., and Terebella lapidaria Linnaeus, 1767 are the main contributors in the dissimilarity of group III with the other groups. Likewise, the dissimilarity-contributing taxa between group IV and the other assemblages were the crustaceans Ampelisca sp Krøyer, 1842, Apseudes sp. Leach, 1814, Caprella acanthifera Leach, 1814, and U. elegans, with N. trispinata and S. armiger (Table 7).Table 7
Results of SIMPER analysis showing the average dissimilarity between the different benthic assemblages. Assemblages were identified by the hierarchical ascendant classification analysis and the species contributing in the dissimilarity of each benthic assemblage (cutoff 50%).
GroupsG I and G IIG I and G IIIG I and G IVG II and G IIIG II and G IVG III and G IVAverage dissimilarity (%)56.8673.2276.0367.0472.7166.01Species contribution, %Abra tenuis3.603.872.78Actinia equina4.233.092.362.60Ampelisca sp.3.022.863.12Amphitoe ramondi2.782.62Ampithoe ferox3.762.782.03Aphelochaeta sp.3.813.662.87Apseudes latreillii1.51Apseudes sp.2.542.422.64Bittium reticulatum2.15Caprella acanthifera3.493.273.57Caprella takeuchii1.651.80Cerastoderma edule3.122.491.73Chone duneri1.911.571.67Cirriformia tentaculata4.384.211.54Clibanarius erythropus2.49Cyathura carinata3.884.782.23Diopatra sp.4.133.233.973.06Dynamene sp.1.53Euclymene palermitana3.133.592.321.713.71Hediste diversicolorHeteromastus filiformis3.705.192.262.40Lysianassa ceratina3.054.651.623.51Lekanesphaera rugicauda3.561.98Leucothoe richiardii2.922.792.17Loripes orbiculatus1.85Lagis koreni2.69Melita palmata2.582.231.641.68Microdeutopus chelifer4.513.623.381.56Nicomache (Loxochona) trispinata2.722.572.80Ophiura sp.2.951.81Peringia ulvae4.468.396.104.983.55Sabella spp.2.431.91Scoloplos armiger2.362.353.06Scrobicularia plana2.803.812.561.84Tanais dulongii2.722.05Terebella lapidaria5.184.952.87Tritia pfeifferi1.68Urothoe elegans3.062.893.16The first two axes of the PCO analyses explained 57% of total variation, and the PCO plot indicated a distinct pattern of the benthic communities’ structure between the “northern” and the “southern” site’s samples (Figure3). They were clearly separated along the PCO1 axis (37.4%) which is negatively correlated with C. tentaculata (r = −0.83), Chone duneri Malmgren, 1867 (r = −0.72), Ampithoe ramondi Audouin, 1826 (r = −0.91), and Ampelisca sp. (r = −0.64) and correlated positively with A. tenuis (0.68), while the PCO2 axis (19.6% of total of variation) separating between groups in both regions is negatively correlated positively with Bittium reticulatum da Costa, 1778 (r = 0.61), Capitella capitata Fabricius, 1780 (r = 0.66), Ampelisca sp. (r = 0.67), Apseudes sp. (r = 0.63), and negatively with Aphelochaeta sp. (r = −0.72).Figure 3
Principal coordinate analysis (PCO) of benthic fauna mean abundance data of sampled sites. Attributes that were most associated to the PCO axes (Person correlations > 0.3) are represented by vectors that overlap the PCO.


## 3.4. Relationship between Environmental and Biological Data

The sequential DistLM analysis showed that the belowground biomass of the seagrass, the salinity, and the percentage of mud fraction in the sediment had a significant correlation on the latitudinal distribution of benthic assemblages (p<0.05), explaining the greatest proportion (63%). However, the best solution provided through the DistLM analysis was found when using five variables (salinity, temperature, percentage of mud, aboveground, and belowground biomasses) as environmental predictors of benthic macrofauna composition, explaining 78% of the total variability between samples (Table 8). The first two dbRDA axes captured 69.8% of the variability in the fitted model and 54.6% of the total variation in the data cloud (Figure 4). The first dbRDA1 axis (36% of the total variation) is correlated with belowground biomass (r = −0.71) and % salinity (r = −0.46). The dbRDA2 axis represents 18.6% of the total variation and correlated strongly with the percentage of mud (r = −0.61) and salinity (r = −0.6) (Figure 4).Table 8
Results of the sequential test of the multivariate regression analysis (DistLM). Environmental variables were analyzed individually (marginal test) and sequentially using a forward selection procedure (adjustedR2 = an information criterion).
VariablePseudo-Fp valueProportionCumulativeproportionBelowground biomass3.83360.00040.32390.3237Salinity2.64180.00830.18520.5092Percentage of mud2.0680.04740.12580.6350Aboveground biomass1.26970.28820.07390.7089Temperature1.36130.26050.07390.7828Figure 4
Two-dimensional distance-based redundancy analysis (dbRDA) ordination representing the model of spatial variation in macrofaunal community structure related to the predictor variables selected through the best linear models based on distance (DistLM).


## 4. Discussion

Latitudinal diversity gradient, peaking in the tropics and declining near the poles, forms the most remarkable large-scale biotic pattern common for both marine and terrestrial systems. Several studies conducted at the regional scale (1000 s of km) have reported a significant change in benthic communities along a latitudinal scale [48, 49]. These changes refer mostly to major effects, such as the proximity of upwellings, variation in water temperature, and the anthropic disturbances [50, 51].Previously, latitudinal comparisons of benthic community structure consider literature reviews and qualitative work [52]. The insufficiency of existing quantitative baselines on this scale limits our ability to assess if changes in these habitats are varying naturally or a result of anthropogenic influence or a mixture of both. Boutoumit et al. [53] found no relationship between latitude and species richness and taxonomic diversity by compiling checklist data considering 12 SECSs along the Moroccan coast which include our studied sites. Our work provides a structured quantitative characterization of Zostera noltei beds in the semi-enclosed benthic ecosystems on the Atlantic coast of Morocco.This study revealed diverse benthic macrofauna for theZostera noltei beds with overall 96 taxa, and crustaceans had the highest abundance compared to the other benthic fauna such as mollusks and annelids. This finding was in line with the results of a study conducted by Tanner [54], which stated that crustaceans were the most abundant group of fauna living in seagrass ecosystems. Comparison of mean species richness with previous studies was difficult because of variance of the sampling design (core dimensions, number of replicates, and number of selected sites). However, given the number of sampling size used in this analysis, the Z. noltei beds of our studied sites (96) were less than those recorded on the same ecosystem in the Kneiss Islands, Tunisia (148), and Arcachon Bay, France (117) [55]. The less diverse Z. noltei beds reported in the different studies where the seagrass is annual [56] or shows large annual fluctuations due to the grazing pressure by migrating seabirds (Table 9).Table 9
List of the taxa found inZostera noltei beds in the study sites.
PhylumClass/orderTaxaCnidariaActinia equina (Linnaeus, 1758)NemerteaNemertea spp.PlatyhelminthesPlatyhelminthes spp. Minot, 1876AnnelidaPolychaetaAphelochaeta sp. Blake, 1991Aphrodita sp. Linnaeus, 1758Capitella capitata (Fabricius, 1780)Chone duneri Malmgren, 1867Cirriformia tentaculata (Montagu, 1808)Diopatra sp. Audouin and Milne Edwards, 1833Diopatra neapolitana Delle Chiaje, 1841Euclymene palermitana (Grube, 1840)Glycera tridactyla Schmarda, 1861Hediste diversicolor (O.F. Müller, 1776)Heteromastus filiformis (Claparède, 1864)Lagis koreni Malmgren, 1866Lumbrineris latreilli Audouin and Milne Edwards, 1833Magelona sp. F. Müller, 1858Marphysa sanguinea (Montagu, 1813)Naineris laevigata (Grube, 1855)Nephtys hombergii Savigny in Lamarck, 1818Nicomache (Loxochona) trispinata Arwidsson, 1906Notomastus latericeus Sars, 1851Owenia fusiformis Delle Chiaje, 1844Phyllodoce sp.Lamarck, 1818Sabella spp. Linnaeus, 1767Scoloplos armiger (Müller, 1776)Spio sp. Fabricius, 1785Streblospio shrubsolii (Buchanan, 1890)Syllides sp. Örsted, 1845Terebella lapidaria Linnaeus, 1767ArthropodaAmphipodaAmpelisca sp. Krøyer, 1842Ampithoe ferox (Chevreux, 1901)Ampithoe ramondi Audouin, 1826Caprella acanthifera Leach, 1814Caprella grandimana Mayer, 1882Caprella takeuchii Guerra-García, Sánchez-Moyano, and García-Gómez, 2001Corophium sp. Latreille, 1806Elasmopus rapax Costa, 1853Ericthonius punctatus (Spence Bate, 1857)Gammarella fucicola (Leach, 1814)Leptocheirus hirsutimanus (Spence Bate, 1862)Leucothoe sp. Leach, 1814Leucothoe richiardii Lessona, 1865Lysianassa ceratina (Walker, 1889)Melita palmata (Montagu, 1804)Metaphoxus simplex (Spence Bate, 1857)Microdeutopus chelifer (Spence Bate, 1862)Urothoe elegans Spence Bate, 1857CumaceaIphinoe serrata Norman, 1867Iphinoe trispinosa (Goodsir, 1843)DecapodaAthanas nitescens (Leach, 1814 (in Leach, 1813–1815))Carcinus maenas (Linnaeus, 1758)Clibanarius erythropus (Latreille, 1818)Hippolyte sp. Leach, 1814 (in Leach, 1813–1815)Pagurus sp. J.C. Fabricius, 1775Palaemon elegans Rathke, 1836Processa sp. Leach, 1815 (in Leach, 1815–1875)Upogebia pusilla (Petagna, 1792)IsopodaCyathura carinata (Krøyer, 1847)Dynamene edwardsi (Lucas, 1849)Dynamene sp. Leach, 1814Idotea chelipes (Pallas, 1766)Lekanesphaera levii (Argano and Ponticelli, 1981)Lekanesphaera rugicauda (Leach, 1814)TanaidaceaApseudes sp. Leach, 1814Apseudopsis latreillii (Milne Edwards, 1828)Chondrochelia savignyi (Kroyer, 1842)Tanais dulongii (Audouin, 1826)DipteraChironomidae larvaeDolichopodidae larvaeMolluscaBivalviaAbra alba (W. Wood, 1802)Abra tenuis (Montagu, 1803)Brachidontes sp. Swainson, 1840Cerastoderma edule (Linnaeus, 1758)Gastrana fragilis (Linnaeus, 1758)Lithophaga lithophaga (Linnaeus, 1758)Loripes orbiculatus Poli, 1795Polititapes aureus (Gmelin, 1791)Ruditapes decussatus (Linnaeus, 1758)Scrobicularia plana (da Costa, 1778)Solen marginatus Pulteney, 1799Venerupis corrugata (Gmelin, 1791)GastropodaBittium reticulatum (da Costa, 1778)Bulla striata Bruguière, 1792Cochlis vittata (Gmelin, 1791)Crepidula porcellana Lamarck, 1801Gibberula miliaria (Linnaeus, 1758)Haminoea navicula (da Costa, 1778)Hydrobia acuta (Draparnaud, 1805)Peringia ulvae (Pennant, 1777)Persicula cingulata (Dillwyn, 1817)Prunum amygdalum (Kiener, 1841)Tritia pfeifferi (Philippi, 1844)Tritia reticulata (Linnaeus, 1758)EchinodermataEchinoideaParacentrotus lividus (Lamarck, 1816)OphiuroideaOphiura sp. Lamarck, 1801The macrofaunal diversity appears to be a more significant structural parameter varying across the sites [57]. It is the most elementary parameter employed by studies examining the large-scale variation of biodiversity and frequently reveals a linear relationship with latitude [58, 59]. Our results state that the species richness of benthic communities associated with Z. noltei was highly variable across the study sites. These results agree with many existing macroecological studies that showed that richness overall decreased with increasing latitude [60].The multivariate analyses of assemblage’s composition showed a 57% of total variation observed in benthic assemblages, with the presence of species that are omnipresent in all sites such asC. edule, C. carinata, I. chelipes, and H. diversicolor. On the contrary, some species were present in just one site (e.g., Lekanesphaera rugicauda Leach, 1814 in Merja Zerga lagoon, Caprella acanthifera Leach, 1814, and Caprella takeuchii Guerra-García, Sánchez-Moyano, and García-Gómez, 2001 in Dakhla Bay). This similar pattern of restrictedness of the species has been stated elsewhere [58]. Furthermore, the PERMANOVA analysis confirmed that faunal assemblages changed significantly among sites.Such variations in benthic macrofaunal composition over a large scale could rise from the supposition that each habitat has its unique characteristics, which recommend an individualistic approach to ecosystem ecology. Certainly, all habitats are vulnerable to environmental and climatic influences, and their variations generate a response from populations [61]. Diversity and spatial heterogeneity of species can be affected by ecological and evolutionary processes at local and regional levels [62], the systemic variability of transitional waters [63], and the resultant distribution of the benthic macrofaunal organisms according to their functional features and niche demands [64].According to DistLM results, variations in belowground biomass, salinity, and percentage of mud explained a larger part of the variation in benthic fauna than other abiotic parameters within the studied systems. Indeed, the presence of seagrass influences the macrofaunal diversity and biomass [65] and the biomass of seagrass affects the organization of benthic macrofaunal assemblages [66]. Moreover, variations might be related to substrate type and organic residues [67].The studied semi-enclosed ecosystems are subjected to different anthropogenic disturbances. Boutahar et al. [68] have showed a clear variation on chemical elements accumulated by Z. noltei leaves along the North-South latitudinal scale of the Atlantic coast of Morocco. Differences in anthropic pressures and environmental conditions do not only exist between sites but also between upstream and downstream stations within the semi-enclosed systems. Furthermore, they are situated along a widespread gradient of climate regimes (semiarid and arid) and hydrological (e.g., temperature and salinity) conditions. The absence of a clear pattern with latitude means that the natural mechanisms that can affect diversity across sites are mostly the same along this stretch of the Moroccan coast. We support the postulate that other different processes, operating at different spatial scales, may explain the latitudinal trends in diversity [9, 69]. Southern systems (Khnifiss and Dakhla) exhibited the highest specific richness and composition variability which can be related to the strong hydrodynamic conditions ensuring the homogenization of water masses and their fast renewal rate. It is also linked with the absence of continental freshwater inputs [35] which maintain the lagoon’s generally good environmental quality [27].


## 5. Conclusions

Our study revealed broad-scale variability in species composition of benthic macrofauna associated toZ. noltei beds along the Atlantic coast of Morocco. The observed variability was influenced by seagrass biomass, which varied greatly crosswise sites. While such patterns could be related to ecological and biological factors such as habitat heterogeneity and traits of component species, regional and biogeographic factors such as climate and oceanographic current regimes may also be important although these variables were not tested directly in this study. The results of our study constitute as baseline data for planning the broad-scale conservation of biodiversity in seagrass beds of Morocco, which remain suffering from multiple human-induced threats such as coastal developments and climate change.



---

*Source: 9199654-2021-10-18.xml*