Latitudinal Variation into the Macrofaunal Assemblages Associated to Zostera noltei Seagrass along the Atlantic Coast of Morocco

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 to Zostera 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.(e ecological community descriptors differ significantly at the level of the site. Specific richness showed an inconsistent significant patternwith latitude.(emultivariate 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. (e present study marked the first attempt on broadscale ecological research of seagrass beds inMorocco 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.


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]. e 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]. e main drivers behind the LDG are focused on theories that are broadly linked to the current climate, historical effects, and biome area [5,6]. e 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]. ey support abundant and generally well-known macrofauna, the density and diversity of which also surpass the existing in nearby bare sediment fields [14,15]. e 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. e dwarf eelgrass Zostera 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]. e 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 with Z. noltei habitats of Moroccan SECSs. en, 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.

Study Sites.
Benthic macrofauna associated to Zostera 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). e 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 km 2 , 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].
e Sidi Moussa lagoon (32°52′N-8°51′W) covers an area of 4.2 km 2 with a maximum depth of approximately 5 m, which decreases progressively towards the upstream part of the lagoon. e 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].
e 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 km 2 . It is characterized by a semidiurnal tide and entrances ranging from 0.8 m to 3.6 [32]. e lagoon is subject to various stressors related to fishing activities, aquaculture, and algae exploitation [33].
e 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 km 2 , and a maximum depth of 8.7 m [34]. e 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].
e Dakhla Bay (23°35′N-15°50′W) extends over 37 km length and 12 km width for a total surface of 400 km 2 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].

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. e 10 replicates per station totalize a surface of 0.12 m 2 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. e hydrological parameters (water temperature, pH, and salinity) were also measured in situ using a HANNA portable multiparameter. For Z. 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%). e 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]. e 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 CO 2 and CaCO 3 percentages after a 1400°C dioxygen burning and a mineral decarbonizing with a sulfuric acid solution [41,42].
Additionally, three replicates per station of Zostera noltei beds were randomly sampled for seagrass characterization: shoot density (nb shoots/m 2 ), aboveground biomass (gDW/ m 2 ), and belowground biomass (gDW/m 2 ). Samples of Z. noltei were cleaned and separated into leaves and belowground parts (roots and rhizomes) and then were ovendried at 60°C until constant dry weight (DW).
Abundance (N: Ind./m 2 ), number of species (S), Shannon-Wiener diversity index (H′, log 2 ), and Pielou's evenness index (J′) were calculated for each sample using e 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]. e 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. e adjusted R 2 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).

Environmental Variables.
Values of environmental parameters are shown in Table 1.
e 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. e sediment composition revealed that the mud content is varying from 8.00 to 88.43% in the sampling sites. Concentration of CO 2 and CaCO 3 showed a significant difference among sites (p < 0.05) with no difference between downstream and upstream samples.
For 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). e 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).

Diversity and Species Composition.
is study identified a total of 17,320 individuals belonging to 96 benthic macroinvertebrate taxa. At the level of the phylum, the Zostera noltei beds in the sampled SECS were highly dominated by crustaceans with 38 species, representing more than 38% of the sampled macrofauna. e 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. e 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 (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). e 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). e other ecological descriptors did not show a significant correlation with latitude.

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). e interaction between effects was significant (p perm � 0.0001) (   (Table 7). e 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 ( Figure 3). ey were clearly separated along the PCO 1 axis (37.4%) which is negatively correlated with C. tentaculata

Relationship between Environmental and Biological Data.
e 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). e 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). e first dbRDA1 axis (36% of the total variation) is correlated with belowground biomass (r � −0.71) and % salinity (r � −0.46). e 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).

Discussion
Latitudinal diversity gradient, peaking in the tropics and declining near the poles, forms the most remarkable largescale 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]. ese 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]. e 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 semienclosed benthic ecosystems on the Atlantic coast of Morocco.
is study revealed diverse benthic macrofauna for the Zostera noltei beds with overall 96 taxa, and crustaceans had the highest abundance compared to the other benthic fauna such as mollusks and annelids. is 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]. e 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).
e 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. ese results agree with many existing macroecological studies that showed that richness overall decreased with increasing latitude [60]. e 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 as C. 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). is similar pattern of restrictedness of the species has been stated elsewhere [58]. Furthermore, the PERMANOVA analysis confirmed that faunal assemblages changed significantly among sites.    Table  4: Results of two-way ANOVA testing for macrobenthic assemblages' differences in the total number of individuals (N), total number of taxa (    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]. e 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. e 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 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%). Groups G I and G II  G I and G III  G I and G IV  G II and G III  G II and G IV  G III and G 20 40 60 Simprof Groups G I G II G III G IV Figure 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). of continental freshwater inputs [35] which maintain the lagoon's generally good environmental quality [27].

Conclusions
Our study revealed broad-scale variability in species composition of benthic macrofauna associated to Z. noltei beds along the Atlantic coast of Morocco. e 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. e results of our study constitute as baseline data for planning the broadscale conservation of biodiversity in seagrass beds of Morocco, which remain suffering from multiple human-induced threats such as coastal developments and climate change.

Data Availability
Data used to support the findings of this work are available, upon request, from the corresponding author.

Conflicts of Interest
e authors declare that there are no conflicts of interest.

Authors' Contributions
Oussama Bououarour, Soilam Boutoumit, and Reda El Kamcha conducted investigation, curated data, wrote the original draft, and performed visualization and Hocein Bazairi conceptualized the study, formulated the methodology, wrote the original draft, reviewed and edited the manuscript, and supervised the work.