Introduction

Nowadays, there is no doubt that reintroduction, defined as attempt to repopulate areas where a species has been extirpated (Kleiman 1989), is accepted as an effective effort in biodiversity conservation (Kleiman 19,890; Seddon et al. 2007; Armstrong and Seddon 2008). Such programs (Moorhouse et al. 2009; Bernardo et al. 2011) are often associated with post-release monitorings to confirm the relevancy of reintroductions efforts. Such post-release monitoring is crucial to: (i) inform on the success of the release operation (IUCN 1998); (ii) know the seasonal pattern of habitat use (Hanane and Magri 2016); and (iii) and (iii) fill gaps in spatial dynamic of reintroduced individuals by recording their spatial position and habitat use.

Since 1942, Morocco has been committed to undertaking in-situ conservation approaches for biodiversity. Therefore, a total of ten national parks, three biosphere reserves, and 154 biological and ecological interest sites (BEIS) were created (Marraha 2007; Ramzi et al. 2018). In this North African country, post-release monitoring programs have been carried-out for some threatened mammal species such as dorcas gazelle (Gazella dorcas Linnaeus, 1788) in M’Sabih Talaa Reserve (Ait Baamrane et al. 2012, 2013), Aoudad (Ammotragus lervia Pallas, 1777) in Amassine reserve (Ramzi et al. 2018), Mhorr gazelle (Nanger dama mhorr Pallas, 1766) in Safia reserve (Abaigar et al. 2019), Atlas deer (Cervus elaphus barbarus Linnaeus, 1758) in Tazekka National Park (Ismaili et al. 2018), and more recently Cuvier’s gazelle (Gazella cuvieri) in Tergou reserve (present study).

According to the IUCN Red List (IUCN 2018), Cuvier’s gazelle is classified as a vulnerable species. In Morocco, this mammal is categorized among protected species by the law on the protection of species of flora and wildlife and the control of their trade (law 29–05, adopted in 2011), and the law on the hunting control (Dahir of July 21, 1923 on the hunting police) (Hanane and Amhaouch, in press). This ungulate is exclusively distributed in three North African countries, including Morocco, Algeria, and Tunisia (Beudels et al. 2013; Gil-Sánchez et al. 2016). The total estimation of the Moroccan natural population ranged between 1900 and 3900 individuals (IUCN 2018; Hanane and Amhaouch, in press). However, during these last years, an attested reduction of the G. cuvieri population is recorded (IUCN 2018), which is mainly attributed to poaching, predation, and habitat loss (IUCN 2018). The reestablishment of this breeding population is thereby necessary with the view to restore its status. As part of its strategic efforts to strengthen wild populations, Morocco’s Water and Forest Department, in collaboration with the Emirates Center of Wildlife Propagation and Al Maha Farm, has carried out a Cuvier’s gazelle reintroduction program in a reserve in the Central-Plateau, the Tergou Forest Reserve (TFR). In this enclosure, a total of 24 Cuvier’s gazelles were released: 12 in November 2018, 6 in December 2018, 3 in February 2019, and 3 in July 2019. Before the release of the gazelles, three management actions were first been undertaken: (i) the establishment of a fence to protect released individuals, (ii) the opening of a track to facilitate their movement, and (iii) the creation of an artificial watering point.

Monitoring of released gazelles in forest protected areas would allow understanding how they react to management actions, especially the opening runway obtained by transforming a portion of the forest habitats (trees, shrubs and herbaceous layer) into sandy open areas. To our knowledge, such assessment has never been undertaken in forest areas (see for exemple, Abaigar et al. 2005 and Moreno et al. 2020). To understand such a pattern, the use of indirect sign sampling techniques is recommended because of their wide use in studying gazelles (Abaigar et al. 2005, 2013; Wronski and Plath 2010; Gil-Sánchez et al. 2016). The number of active latrines (NAL) are among the indirect signs that are highly used (Wronski and Plath 2010; Herrera-Sánchez et al. 2020). It is a relevant tool for surveying gazelle populations as has been acknowledged by Wronski and Plath (2010). According to Gosling and Roberts (2001) and Wronski and Plath (2010), latrine surveys can also inform on (i) territorial organization, (ii) social behaviour (Wronski et al. 2010), and (iii) population size (Wronski et al. 2013). Their use is also favored by easy recognition (Wronski and Plath 2010; Wronski et al. 2013) and good availability.

The assessment of spatio-temporal latrine distribution provide information about distribution and abundance of populations, and are used as the basis for management decisions (e.g., Patterson et al. 2004; Crowley et al. 2013). Indeed, knowledge of areas with strong and weak densities of latrines remains a relevant tool promoting the decision making by managers for a better conserve reintroduced gazelle populations. By knowing spatio-temporal patterns governing latrine distribution, the establishment of new populations of Cuvier’s gazelle populations will be possible under better auspices.

In this study, we aimed to evaluate the effects of season, microhabitat and management strategies (i.e. artificial watering point, track, and fence) on the spatio-temporal variations of NAL in the TFR. The results of this investigation might be useful to improve the management and conservation of this vulnerable mammal species.

Materials and methods

Study area

The study was conducted in the TFR, a 118-ha enclosed forest located in the Central-plateau, in north-central Morocco (33°26′34” N, 6°7′39″ E). The TFR, located in a steep mountain side between 544 and 1053 m a.s.l., was created in 1998 to conserve the ungulates that were historically present in north-central Morocco. The study area has a subhumid climate characterized by wet and relatively cold and rainy winters (average = 773.2 ± 25.1) and hot and dry summers (average = 38 °C ± 2.11) (Cuzin 2015b; Haraz 2020).

A water tower was established inside the TFR to ensure the availability of water for the reintroduced gazelles. A runway was also managed to connect the summit and the base of the mountain side. The fence surrounding the TFR is 5024.76 m long. In this reserve, the forest stand is dominated by Olea europeae, Quercus ilex, Tetraclinis articulata, Buxus balearica, Phyllirea latifolia, Cistus salviifolius, Cistus creticus, and Chamaerops humilis. Around the TFR, the landscape is composed by the (i) same forest type, (ii) cerealiculture especially Barley (Hordeum vulgare), and (iii) culture of an aromatic and medicinal plant, i.e., lavender (Lavandula angustifolia).

The TFR is relatively isolated from the humans, with the nearest habitation from the fence located at a distance of 3.25 km.

Data collection

Field surveys

The data on the NAL were collected monthly from June 2019 to February 2020 along six 500-m transects (it was not possible to collect data between March and June 2020 because of the lockdown established in Morocco du to COVID-19 pandemic). This number of transects was chosen according to the availability of human and material resources. The total transect length was 3 km and the area sampled was 3 ha (~3% of the total area of the TGR); all latrines present within 5 m on each side were considered. As transect lines should be placed at random with respect to the distribution of the animals (Marshall et al. 2008), six points, corresponding to starting points for each transect, were randomly selected using the Quantum GIS (QGIS) random selection tool (Quantum GIS Development Team 2017) (Fig. 1). Then, they were fixed to count the NAL monthly. We collected the NAL (active = presence of fresh fecal deposits; inactive = absence of fresh fecal deposits) per transect at each monthly passage along the six transects. This was relatively easy due to the sole presence of Cuvier’s gazelle in this enclosure. The direction of the passages was always north-south with no return. Coordinates of all observed NAL were recorded using a global positioning system (GPS) device (Garmin eTreex HC) (± 3 m).

Fig. 1
figure 1

Map showing the geographic location of the Tergou reserve in Morocco and the position of the six transects in this protected area (Hanane and Amhaouch 2021)

Monthly observations were pooled and presented as seasonal observations as follows: summer (end June to August), autumn (end September to November), and winter (end December to February). Observation periods were restricted to daytime hours in the three seasons: 06:00–20:00 h (summer), 07:00–19:00 h (autumn), and 8:00-18:00 h (winter).

Predictor variables

A total of eight variables related to season, microhabitat features and management actions were selected (Table 1). These variables were used because they can potentially influence the NAL in this protected area. Indeed, the number of gazelles’ latrines is known to varied according vegetation (Attum et al. 2006; Attum 2007; Wronski and Plath 2010; Attum and Mahmoud 2012), season (Abaigar et al. 2005, 2013) and management variables, mainly distance to the fence (Abaigar et al. 2005, 2013). All distance variables (i.e., management variable (Table 1)) were measured using QGIS, whereas microhabitat features were collected in circular plots of 10-m radius (0.03 ha) (Maziarz and Broughton 2015). Within each of these plots, we visually estimated the covers (%) of trees, shrubs, herbaceous layer, and bare ground (Table 1) (Alaya-Ltifi and Selmi 2014; Hanane et al. 2019).

Table 1 Description of variables with their means, standard deviations, and ranges (Min. –Max.). Microhabitat variables were collected within circles of 10-m radius. Variance Inflation Factor (VIF) for each variable is added

Statistical analysis

Statistical analyses were performed in R v4.0.3 (R Core Development Team 2018). To ensure normality, one independent variable (distance to the nearest fence) was log-transformed. To overcome multicollinearity, all predictive variables were tested with the Variance Inflation Factor (VIF) (Quinn and Keough 2002) using the R package “car” (Fox and Weisberg 2011). Variables with a VIF > 3 were removed as recommended by Zuur et al. (2010). To do this, we gradually removed the variable with the highest VIF until the set of variables had a VIF lower than 3 (Table 1).

To test the effects of season, microhabitat features and management variables on NAL, Generalized Linear Mixed Models (GLMM) with a Poisson error distribution was performed using the package “lme4” (Bates et al. 2014). Transects were included as random factor in the model to account for potential non-independence of multiple observations at the same transect. An all-inclusive design (all possible combination models) was developed using Information-Theoretic approach (Burnham and Anderson 2002). Models were then ordered by increasing Akaike Information Criterion corrected for small sample sizes (AICc) (Burnham and Anderson 2002) using the package “MuMIn” (Bartoń 2015). All models with ΔAICc lower than 2 were considered as equally good (Burnham and Anderson 2002). Variance explained was calculated using the method of Nakagawa and Schielzeth (2013). Marginal and conditional R2, respectively, describe the variance explained by fixed effects and the full model, were calculated using the function “rsquared.glmm” in the R package “MuMIn”.

To ensure that observations were independent and to be able to address spatial autocorrelation in data, we used the indicator semivariogram (Goovaerts 1998) on the residual of the best model in terms of AICc value. 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 is considered as strong spatial dependence, while 0.25 < NSR < 0.75 and NSR > 0.75 indicate moderate and weak spatial dependence, respectively (Cambardella et al. 1994). The packages “sp” (Pebesma and Bivand 2005), “lattice” (Sarkar 2015) and “gstat” (Pebesma 2006) were used to draw semivariograms. When spatial autocorrelation was encountered, we used spatial generalized linear mixed models, fitted via penalized quasi-likelihood (glmmPQL) using the package “MASS” (Venables and Ripley 2002), which enables the construction of spatial models with dependent data that are not normally distributed and which is among the best techniques for this kind of data (Dormann 2007). We adopted an exponential spatial correlation structure; however, tests with gaussian and spherical structures led to the same results. We also created spatial variables by means of Moran’s eigenvector maps (MEMs), using the packages “spdep”, “adespatial”, “ade4”, “maptools” and “adegraphics” (Dray et al. 2006), a method that produces flexible spatial predictors through principal coordinate analysis of a truncated geographic distance matrix among different points, while capturing spatial effects at multiple spatial scales. These orthogonal spatial variables (MEMs) can be used in regression (independent variables) to account for spatial autocorrelation (Assandri et al. 2019; Farhi et al. 2020) and for variation partitioning to account for the space (Hanane 2018; Farhi et al. 2020; Hmamouchi et al. 2020).

For subsequent variation partitioning (VP) analyses, we retained only the variables composing the final best model AICc with confidence intervals of parameter estimates not encompassing zero. VP was applied to evaluate the specific contribution of season, microhabitat features and management variables and their joint fractions in the NAL. VP was conducted by including, in a first step, all final variables of each group constituting the best model AICc: (i) season, (ii) microhabitat, (iii) management and (iv) space (MEM), and in a second step, only three classes of variables: (i) the interaction term between season and management variables (see results), (ii) microhabitat, and (iii) space (MEM). Consideration of the interaction effect was necessary to quantify its contribution for explaining the NAL. We tested for the significance of the unique fractions using the function “rda” from the vegan package (Oksanen et al. 2013). However, it was not possible to test the significance of the shared variation (Truchy et al. 2019).

We used the package ‘visreg’ (Breheny and Burchett 2012) to plot the relationship between the number of active latrines and the interaction between season and distance to the nearest developed runway.

Results

In the Tergou reserve, the maximum number of active latrines (NAL) per transect was four in summer (2019), six in autumn (2019), and three in winter (2020). The average characteristics of those latrines are summarized in Table 1.

The top-ranking model assessing factors influencing number of active latrines indicated that there is significant variation of this number according the interaction between season and distance to the nearest developed runway. Although included in the best AICc combinaison (Table 2), the distance to the nearest fence is not very relevant in explaining the NAL (Table 3). Nugget-to-sill ratio (NSR) showed that the residuals of these best models AICc exhibited strong spatial autocorrelation (Table 2). When we considered spatial autocorrelation in the modeling through glmmPQL models (Table 3), the effect of interaction between season and distance to the nearest developed runway remained significant, whereas that of the distance to the nearest fence became marginally significant (Table 3).

Table 2 Competitive models explaining active numbers of latrines per transect
Table 3 Non-spatial and spatial model parameters for the two models explaining the number of active latrines per transect in the Tergou reserve, Morocco. See methods for details

In the enclosure of Tergou, the number of active latrines was affected by the interaction between season and distance to the nearest developed runway. Indeed, glmmPQL model showed that in autumn, the NAL was greater at the farthest distances from the developed runway (Fig. 2), whereas, in the summer and the winter, this number was rather greater in proximity to the developed runway (Fig. 2). This same spatial model also showed a marginally contribution of the distance to nearest fence in explaining the NAL (Table 3), being closer to it whatever the season.

Fig. 2
figure 2

Relationship between the number of active latrines (NAL) per transect (predicted values) and the distance to the nearest developed runway in the Tergou reserve (Morocco) during summer (dark gray line), autumn (light gray line) and winter (black line)

When we considered each set of variables alone (i.e., season, microhabitat, management and space (MEM)) in explaining the NAL, the VP analysis revealed that the pure effect of season was robust (Adj. R2 = 0.31, p = 0.001) (Fig. 2a), followed by the joint effect of season and management (Fig. 3a) and the unique fraction of management (Adj. R2 = 0.09, p = 0.01). In order to quantify the share contribution of the significant effect of interaction between the season and management variables in the variation of the NAL, we considered the terms of this interaction (i.e. season × distance to the nearest developed runway) as a single entity in the VP analysis, in addition to those of microhabitat and space (MEM). This analysis indicated that the interaction effect alone explained 82% (p = 0.001) of the total variation in the NAL (Fig. 3b), followed by the joint effect of the “interaction term” and space (Fig. 3b).

Fig. 3
figure 3

a Venn diagrams for variation partitioning showing the percentage contribution of season, microhabitat, management, and spatial components (MEM) in explaining the NAL. b Venn diagrams for variation partitioning showing the percentage contribution of season × distance to the nearest developed runway, microhabitat, and spatial components (MEM) in explaining the NAL

Discussion

In the present study, we investigated the effects of season, microhabitat features, and management strategies on the spatio-temporal patterns of NAL for the reintroduced Cuvier’s gazelles in the Tergou reserve. Our results revealed that this number depends on the interaction between season and distance to the nearest developed runway, being more important close to this latter during summer and winter and far away from it during autumn. A marginal effect of distance to the nearest fence was also recorded. These results demonstrate the presence of a temporal pattern in the distribution of reintroduced gazelles in the TFR, depending on the distance to the developed runway, and to a very less degree to distance to the nearest fence in every season.

Statistical analyses show that high NAL were recorded on or in the near vicinity of developed runway in summer and winter. This developed runway made possible the availability of sandy open areas on a portion of the hill of the TFR. Such habitat is often desired by some gazelle species like the dorcas gazelles (Abaigar et al. 2013) and mountain gazelles (Gazella gazella Pallas, 1766) (Wronski and Plath 2010), which use, among others, sandy grounds for the latrines deposits. The selection of this habitat can potentially be explained by the need (i) to ensure a high visibility as already evoked by Abaigar et al. (2005), which is facilitated by the presence of a relatively strong hill, and (ii) to have a relatively aerated areas, which relieve the gazelles in periods of strong heats.

The spatial pattern of NAL showed a significant change in autumn. Indeed, contrary to summer and winter, gazelles much more use the latrines located far away from the developed runway. Such a behavior would be related to the availability of sandy ground in the Tergou enclosure. Indeed, in winter, the developed runway is the only area where the sandy ground is available. The winter precipitations do not allow keeping the sandy ground present outside the runway. In summer, although this soil type is available beyond the developed runway, the strong heats generally force the reintroduced Cuvier’s gazelles to use latrines located on or in close vicinity of the track (Kacem et al. 1994; Abaigar et al. 2005). In contrast, in autumn, the use of latrines far away from the track is made possible because of significant decrease of temperature. Our findings were consistent with Abaigar et al. (2013) who recorded a significant high number of latrines in the plateaus characterized by more than 60% of sandy soil for dorcas gazelles in Senegal.

In the Tergou enclosure, microhabitat variables are not robust enough to explain the NAL (presence in the fourth best model AICc, see Table 2). It is possible that because of the habitat homogeneity in this enclosure, no microhabitat variables has emerged as a relevant predictor of the NAL. Overall, in the TFR, Cuvier’s gazelles do not prefer install their latrines in areas with high density of trees and shrubs as it has been also recorded in Tunisia (Abaigar et al. 2005).

Because of its marginally significance, the distance to the nearest fence could be considered as potential predictor explaining the NAL in the TFR. In this enclosure, this number has a trend to increase with proximity to the fence. This is in agreement with several studies which have reported the presence of high abundance of latrines along the boundary of their home range/territory (e.g., Brashares and Arcese 1999; Stewart et al. 2001; Wronski and Plath 2010). The same pattern in latrine distribution was recorded in Senegal for another reintroduced gazelles species, i.e. the dorcas gazelle (Abaigar et al. 2013). The recorded trend in Cuvier’s gazelles regarding the fence would be related to their birthplaces. Indeed, gazelles born in captivity generally consider the fence as a territorial reference (Abaigar et al. 2005), resulting an attraction towards it.

When considering each group of variables alone of the best model AICc, variation partitioning analyses evidenced the robustness of the pure effect of season in explaining the NAL. The shared effect of season and management (15%), and that of the season, management and space (12%) contribute also but with a lesser degrees. In contrast, after including the interaction between season and management as a single entity in the VP, the pure effect of this interaction turns out to be very robust (82%). This reinforce the fact that the effect of distance to developed runway on the NAL is strongly depending on season. This is all the more true that only 8% of the variance remained without no explanation when we have considered the interaction term as a single entity in the VP analysis (versus 35% when considering each group of variables alone).

Conclusions, implications and perspectives

This study, conducted in a Central plateau forest reserve, allowed a better understanding the distributional dynamic of Cuvier’s gazelles using the NAL as dependant variable. Our results suggest that the NAL is robustly dependent on the interaction between season and management variables (i.e., distance to the nearest developed runway), and, to a lesser extent, of distance to the nearest fence whaterver the season.

For the species’ management and conservation, the presence of sandy grounds in Cuvier’s gazelles reintroduction areas is probably to be considered in the view of a release in the wild in this forest area of Moroccan Central plateau. However, before preconizing such a release in the wild, it would likely be relevant to extend the surface area of the TGR further as well as those of sandy grounds. This would guarantee ideal conditions for gazelles particularly in summer and autumn. Such need would, however, useless in arid protected areas where sandy grounds are well represented.

Furthermore, it is recommended to test the geographical generality of our results using the same analytical approach in other protected areas and in the wild to find out if the NAL follows a similar spatio-temporal pattern to that found in the TGR. This would undoubtedly guide future conservation actions for the species both in forest reserves and unprotected forests.

As any enclosure, the TFR would offer a good opportunity for consolidating and improving the knowledge on this vulnerable mammal species in the view of a successful release in the wild. Because the scientific monitoring is crucial for wildlife management (Wronski et al. 2013), post-release research programs must be continued in this enclosure, expecially after its possible extend, with a view to better understand ecological mechanisms governing the distributional dynamics of Cuvier’s gazelles. An attention should also be given to study (i) seasonal diet variation, (ii) composition and structure of vegetation, and (iii) gazelle health (Ezenwa 2004; Abaigar et al. 2013).