Geographical, Environmental and Anthropogenic Factors Affecting the Occurrence of Wintering Eurasian Spoonbills (Platalea leucorodia) in Morocco

Investigating habitat selection and ecological factors trade-offs is a major avian ecology topic which is closely implicated for conservation purposes. Studies dealing with the impacts of ecological factors on wintering Spoonbills are overall scarce. Here, we used Principal Component Analyses (PCA) and Generalized Linear Mixed Models (GLMM) to test the relevance of geographical, environmental and anthropogenic factors in predicting the occupancy and abundance of the Eurasian spoonbills (Platalea leucorodia) during 2008-2011 within 28 Moroccan wetlands. The Eurasian Spoonbill mean annual occurrence was 59.2% (± 1.5% SE; 95%CI: 56.3%-62.1%). Among the occupied wetlands, 15 (83.3%) were regularly occupied. We found that the Eurasian spoonbill occurrence was negatively related to distance to coastline, altitude and human presence, whereas its abundance increased significantly with increasing mudflat areas. These findings highlight the significant effect of human presence in selecting wintering sites, but not in the prediction of abundance. Management strategies should therefore address specific attributes of coastal wetlands and should carefully consider the effects of habitat use especially those related to mudflats. We further suggest guidelines for future studies to understand the dynamic of Eurasian spoonbills wintering in the region.


Introduction
When choosing a site at which to settle, individual waterbirds need to simultaneously consider the quality of available sites in terms of their environmental conditions, resource availability and the presence/frequency of predators and/or parasites [1,2,3,4,5]. However, this choice can be also influenced by: (i) geographical and topographical factors, which may either decrease or enhance attractiveness of a site [6, 7, 8, [9], and (ii) human activities, which are known to affect wildlife populations ( . Patterns of avian habitat use have recently generated numerous interesting questions in community ecology [12][13], and the assessment of avian habitat relationships has become an important part of wildlife, resource management and conservation planning [14,15,16], especially within the context of global change [17]. Wetlands are habitats of the highest conservation interest due to their high biological diversity [18,19,20]). They also play a vital role in the life cycle of many species, both during wintering, breeding and migration of waterbirds [21]. Among them, the Eurasian Spoonbill Platalea leucorodia (hereafter spoonbill) is one of the most relevant species in the Palearctic. It is a widespread wading bird that breeds patchily from Europe to East Asia, India, the Red Sea and Northern Africa [22]. The species is labelled to be of 'Least Concern' globally according to IUCN criteria. The global population is of about 66,000-140,000 individuals [23].
The Kingdom of Morocco, located in the extreme north-west of Africa, is bordered to the north by the Mediterranean Sea, to the west by the Atlantic Ocean, and to the south by Mauritania. It supports habitats ranging from high-altitude moorland through cork-oak forests to wetlands, deltas, arid steppes and deserts [24]. Morocco is located at the crossroads of several avian migratory routes between Europe and Africa. Morocco plays a major role in the migration of waterbirds as it supports numerous wetlands and 70% of its 3500 kilometers of coastline lie within the East-Atlantic Flyway [25].
In Morocco, the spoonbill is a common passage migrant, winter visitor and localized breeder in Smir lagoon and Oued EL Maleh [26][27]. In the winter period, the species is mostly found in wetlands along the Atlantic coast, especially at estuaries and lower marshes [28]. The Spoonbill is a protected species according to 29.05 national law.
Most studies on Spoonbils have focused on migration [29,30,31,32,33,34], diet [35,36,37,38,39], breeding [40,41,42,43,44], and molecular biology [45][46]. To date there are few studies on abundance and occupancy of wetlands during the winter period in Morocco and even worldwide [47][48]. Yet, there is a real need to understand the ecological factors and processes affecting these two parameters. Such investigations are fundamental to our understanding of the dynamics of spoonbill communities in North African wetlands and in devising conservation plans. Thus, the major objectives of our study were (i) to identify the most important factors affecting wetland occupancy and abundance of spoonbills during winter period (ii) to disentangle the relevance of environmental versus anthropogenic factors by considering wetlands and habitat features along with human presence. Given the high sensitivity of the Eurasian spoonbill to habitat conditions and its low sensitivity to human disturbance [49], we hypothesized that the distribution and abundance of this species in Morocco are largely dependent on wetland and habitat features rather than human presence.

Methods
The 28 wetlands of international conservation importance were studied our of 38 based on their designation as both Ramsar and Important Bird Areas (IBA) sites (www.birdlife.org; n = 20) (Ramsar2021; Ramsar Convention Bureau 2002, Fig. 1) or as exclusively IBA sites ( [50]; n = 6; Fig.1). Among these wetlands, four are of national importance since they are designated, by the Moroccan Protected Areas Study [51], as Sites with Biological and Ecological Interest (SBEI). Eighteen wetlands are coastal areas, while twelve are located in inland areas (Fig. 1).

Occurrence data
Eurasian Spoonbill surveys were conducted during four consecutive years, from 2008 to 2011 in December-January because variable used in this study were only collected in this period during the Eurasian Spoonbill censuses, which is in line with the International Waterbird Census (IWC; [52]. All study sites were visited once a year because of limited human resources, shortness of the survey period and the large distances between studied wetlands (up to 1533 km). During each visit, direct counting of spoonbills was carried out using binoculars and a spotting scope from the shore and from accessible vantage points. The number of counting points depended on wetland shape and area size. The number of observation points per site ranged from 3 to 6. Observations were made between sunrise and midday, and in good weather, ensuring that counting conditions were relatively homogeneous.
A wetland was considered to be occupied when at least one individual of Eurasian Spoonbills was recorded. All data were collected by the same author [I. Cherkaoui (IC)].

Geographical, environmental and human data
For each studied wetland, we collected data related to geographical characteristics, habitat features, and human presence ( Table 1). The visual inspection of the resulting proportions revealed no difference between years; data of 2011 were therefore chosen to represent the proportional coverage of habitats. The choice of variables was motivated from our overall aim, which was to identify factors driving wetland site occupancy and abundance of wintering spoonbills in Morocco. Some variables have been already correlated to the occurrence of coastal wintering spoonbills [49].
The location of each wetland (latitude and longitude) and its altitude were determined using a Global Positioning System (GPS) device (Garmin eTrex HC) and compared to altitude data from Google Earth to check their accuracy. We took the wetland surface area from [25] for Ramsar sites, from Birdlife datasets (www.birdlife.com) for Important Bird Areas (IBA), and from the Moroccan Protected Areas Study [51] for SBEI.
To estimate the proportion areas covered by the shallow water and mudflat, we used Google Earth Pro (2017 cover), since high resolution images are available for all studied wetlands. The analysis process was realized as follow: (i) for each wetland, and to bypass any difficulties in identification, many available images were explored to draw polygons for all habitat types; (ii)we proceed in the same way to draw lines for low or high tide; (iii) once finished (for all wetlands), we exported all these polygons and lines to Quantum GIS v2.10; (iv) in QGIS, for each wetland, we gathered, for each habitat type, all its polygons in dedicated layer with specific color; and created the intertidal area polygon by taking into account the lowest and the highest exported lines for tides; and (v)the proportion of the area covered by the shallow water and mudflat, at low tide, were then calculated. In addition, the vegetation cover was calculated from the Google Earth Pro and QGIS according the formula: the ratio of vegetation coverage area to the wetland area × 100.
The distance to the nearest coastline was determined in QGIS using the centroid of each studied wetland. To assess the impact of human disturbances, we collected, for each study year and at each wetland, the number of people present. This has allowed us to determine the human density by dividing the number of persons per surface-area unit.

Statistical analysis 3.1 Preliminary analyses
Statistical models were implemented in R (version 3.1.0); [53]. Before performing statistical analyses, we checked for normality and homogeneity of variance for all variables. In order to approach residuals from normality, variables that did not conform to the requirements for parametric tests were log or square-root transformed prior to all analyses [54,55,56]. We also checked for possible correlations among variables by using Pearson's correlation (r) index. Based on two Principal Component Analyses (PCA) [one for data of occupied and unoccupied wetlands (PCAoc), and the otherone for data of occupied wetlands only (PCAabto model the abundance)], we reduced the total number of studied variables (n = 9) to a smaller number by merging different categories, while minimizing the loss of independent information. We used this reduced set of components (four for PCAoc and four for PCAab) as final variables in our further analyses. The Kaiser-Meyer-Olkin measure of sampling adequacy (KMO) indicated that our data were suitable for the two PCAs (PCAoc: n = 9, KMO = 0.47; Bartlett-test for sphericity, χ2 = 749.18, p< 0.001; PCAab: n = 9, KMO = 0.59; Bartlett-test for sphericity, χ2 = 315.35, p< 0.001, [57].

Model development
Two types of models were used to analyze the relationships between the four PCA components and wetland occupancy and abundance of wintering spoonbills. First, we applied generalized linear mixed models (GLMMs) with binomial error distribution and a logit link function using the lme4 package (version 1.1-7; [58] to test the relationships between wetland occupancy (response variable) and the four explicative PCAoccomponents (independent variables).Second, to determine whether our PCAab components are related to abundance of spoonbills, we used another Generalized Linear Mixed Models (GLMMs) with a negative binomial error distribution and log link function using the lme4 and MASS packages [59]. In both sets of models, we included year and site as random effects to statistically control for betweenyear variation and potential non independence of multiple observations at the same site.
We limited the number of variables included in our models due to the small sample size of wetland sites. Models were then ordered by increasing Akaike Information Criterion (AIC) and AIC weights were calculated [60]. Model with the lowest AIC was selected as the best fitting model. We corrected AIC for small sample size using AICc [60] using theMuMIn package [61]. The top ranked models (within two points of AICc) were then averaged as implemented in the MuMIn R-package [62]. Model weights were used to define the relative importance of each explanatory component across the full set of models evaluated by summing weight values of all models that included the explanatory component of interest.

Spatial autocorrelation
One of the assumptions of parametric statistics is that observations are independent of each other. This assumption is often violated with spatial data, so it is important, when present, to test for and subsequently address spatial autocorrelation in data prior to data analysis. The presence of spatial autocorrelation was tested using semi-variograms of the residuals of the best model in terms of AICc value using the "sp", "gstat", and "lattice" packages [63]. Nugget to total Sill Ratio (NSR) was expressed as the percentage of total semivariance and was used to define for spatial dependency: NSR < 0.25 indicated strong spatial dependence, 0.25 < NSR < 0.75 indicated moderate spatial dependence, and NSR > 0.75 indicated weak spatial dependence [64].
Strong correlations were found between some of the spatial, wetland, habitat and human variables. The PCA summarized these studied field variables into four independent axes accounting for 82.5 % of the variance of the original dataset [28.6 % (eigenvalue = 2.576),22.2 %(eigenvalue = 1.997), 18.2% (eigenvalue = 1.638), and 13.5% (eigenvalue = 1.211) respectively].The varimax rotation revealed a first gradient (PCwetland features) characterized by high loadings of variables related to wetland features as it was positively correlated with distance to the nearest coastline (r= 0.983,p< 0.001), and altitude (r= 0.964,p< 0.001). The second gradient (PCgeographic position) represents an axis of the spatial positionof wetlands as it was positively correlated with longitude (r= 0.955,p < 0.001), and latitude (r=0.950,p< 0.001). The third axis represents an axis of increasing habitat features as it was positively correlated with the shallow water cover (r = 0.801, p< 0.001). The fourth axis represents an axis of increasing human presence as it was positively correlated with human density (r = 0.864, p< 0.001).
The GLMM testing for effects of the four principal components showed that the two models were equally good (  Overall, the wintering spoonbills have especially occupied wetlands located in close proximity to the coastline, at low altitude, and having a low human presence (Fig. 2 and3). With regard to spatial autocorrelation, semi variograms, established from the residuals of the best models, indicate the absence of a spatial pattern in the data (the NSR ratios were equal to 1 for the two top models).

Abundance
In evaluating the effect of ecological factors on the abundance of spoonbills, we used the 64 cases of presence as a dependent variable. As for wetland occupancy, strong correlations were found between some of our studied variables. The PCA summarized these field variables into four independent axes accounting for 87.5 % of the variance of the original dataset [31.5 % (eigenvalue = 2.523),28.3 %(eigenvalue = 2.566), 18.2% (eigenvalue = 1.638), and 13.5% (eigenvalue = 1.211) respectively].The varimax rotation revealed a first gradient (PCwetland features) characterized by high loadings of variables related to wetland features as it was positively correlated with distance to coast (r= 0.923,p< 0.001), altitude (r= 0.910,p< 0.001), and wetland area (r = 0.703, p< 0.001). The second gradient (PCgeographic position) represents an axis of the spatial position of wetlands since it was positively correlated with longitude (r= 0.968, p<0.001), and latitude (r=0.945, p< 0.001). The third axis represents an axis of increasing habitat features since it was positively correlated with the mudflat cover (r = 0.956, p< 0.001). The fourth axis represents an axis of increasing human presence since it was positively correlated with human density (r= 0.972, p < 0.001).
The GLMM testing for effects of the four principal components indicated that the two models were equally acceptable ( Table 2). The most important predictor explaining the abundance of wintering spoonbills was PChabitat features (β= 0.78± 0.18, 95% CI: 0.427 to 1.333; importance =0.994). The PChuman presence (β= 0.13 ± 0.15, 95 % CI: 0.164 to 0.424, importance =0.427), although present in the top model, its contribution is not significantly important [∆AICc = 0.40 (< 2)] ( Table 2).Overall, wintering spoonbills occurred more abundantly in wetlands with large covers of mudflats (Fig. 4). Regarding the spatial autocorrelation, the mean NSR ratio was 0.754 ± 0.06 (95% CI: 0.64 to 0.87). This mean value is situated at the lower limit of the third class established by [64] suggesting a weak spatial pattern in the data.

Discussion
The aim of this study was to assess the relevance of wetland and habitat features and human presence as predictors of the abundance and occupancy of wetlands by spoonbills wintering in Morocco.
Results of the modelling showed that both wetland features and human presence negatively affect the probability of presence of wintering spoonbills in Morocco. In fact, wetlands in proximity to coasts and located at low altitudes are the most occupied by this wading species. This is not surprising to the extent that, in Morocco, this species is mostly found, during the winter period, in wetlands situated on or near coasts especially at estuaries and at lower marshes [28]. This winter distribution could be dictated by the migratory flyway used by this species: along the Atlantic coast between France and Senegal [65]. Indeed, our results support a scenario where wintering spoonbills are reaching suitable areas in Morocco close to the autumnal migratory flyway. It turns out that these suitable areas are in majority coastal Ramsar sites whose importance for waterbirds is no longer to be confirmed . On the other hand, more than 80% of study sites were regularly occupied by wintering spoonbills. Although the present study was carried out over a relatively short period (four years), this rate highlights an important use of these wetlands during the wintering period. Yet, we ignore if the wintering spoonbill are showing fidelity to studied wetlands as we have not tagged birds. But, it seems that this is the case considering this important rate of use. [33] have shown that adult birds are highly site faithful during winter.
The occurrence probability of spoonbills in Moroccan wetlands is also affected, but with a lesser extent, by human presence. This is all the truer that this biotic factor is known to affect the distribution of waterbirds ( [67,68,49]. Indeed, in recent years, increasing human pressures and activities along coastal areas were at the origin of reducing the quality of intertidal habitats in many regions of the world [69][70]. In Morocco, the occupancy of wetlands by wintering spoonbills seems to be a trade-off between the proximity to coastlines where the migration happens, and the human presence. Regarding the abundance, statistical analyses showed that the main factor influencing the number of spoonbills wintering in Morocco is the coverage of mudflats. This is consistent with [66], who demonstrated that in general waterbirds abundance was significantly related to habitat characteristics. [22] have also reported the preference of spoonbills to this type of habitat. In a recent study, conducted in the Tunisian Gulf of Gabès [49], have shown that habitat features, especially large mudflats, have a positive effect on the abundance of Eurasian spoonbills. This trend seems to be largely due to higher food abundance [71]. The effect of large mudflats on the abundance is not limited to the spoonbill, but concerns also other groups of waterbirds like shorebirds and wading birds (e.g. [72,73,74,75]. Surprisingly, human presence had no impact on abundance of wintering spoonbills leading to a key question: why the human presence impacted the probability of presence of spoonbills wintering in Moroccan wetlands and not abundance? We suspect that, once a wetland is occupied, the level of human presence remains tolerable by the spoonbills. This is most likely why the human presence has not emerged as a predictor in the abundance analysis. Another possible explanation is that avoidance of human presence is secondary to selection of habitat. This has been also reported by [76], who have stated the importance of habitat selection in comparison with human presence. Many waterbirds are willing to tolerate increased levels of human disturbance in response to the favorable habitat found [76]. This was confirmed by a study in the Gulf of Gabès (Tunisia) where the abundance of wintering spoonbills was also affected by the habitat features and not by human presence [49]. In Morocco, like in Tunisia [49], the spoonbills occurred mainly in large mudflats, which are in general located far away from humans.
Overall, it appears that, for wintering spoonbills, there is a first level of selection, which bears on the choice of wetlands. This latter takes into account both the wetlands position relative to coasts and human presence. Once settled, the number of spoonbills mostly depends on habitat characteristics such as the proportion of mudflats, while tolerating some level of human presence.

Conclusions and future prospects
The results of this study increase our understanding of the factors shaping the occupancy and abundance of spoonbills wintering in Morocco. Coastal wetlands most likely represent key areas for wintering spoonbills in this region. At the meantime, a low human presence is the guarantor of the presence of the species. Large mudflats play also a key role in increasing their abundance. Occupied wetlands should receive special priority and protection efforts in preservation programs. It is also crucial to maintain a low level of human presence within mudflats. Indeed, considering the Moroccan protected area context, it is important to integrate, in the management plan and the wetland zoning, mudflats as fully-protected habitats given their importance and their usefulness for feeding wintering waterbirds. In parallel, it is of a primary importance to raise awareness among users of these wetlands on the interest of mudflats, not only for the spoonbill but also for other mudflat-dependent species, such as wading and wader species.
Furthermore, additional detailed investigations are needed to allow deeper insights into the presence of spoonbills during winter in Morocco. Our future challenge is to determine quantitatively, via appropriate statistical analyses, the tolerance threshold of human presence under which spoonbills occupy a wetland, but also the tolerance threshold from which the abundance of spoonbills is affected. These recommendations cannot be implemented without the continuation of counting birds while collecting data on human activities, habitat characteristics, food supply, water quality (limnological data), including the conservation status of sites. The Moroccan wetlands, in majority located at the crossroads of several avian migratory routes between Europe and Africa, offer a great opportunity to document the population dynamics of spoonbills both during the winter period and in the course of spring and automn migration.