Introduction 

The Mediterranean biome, which houses 20% of the world’s total floristic richness, has been designated by the IUCN as a global biodiversity hotspot ([22]). Within the Mediterranean basin, the Baetic-Rifian complex, comprising Andalucia (southern Spain) and the Rif mountains (northern Morocco), was recognized as an important hotspot of plant diversity ([39]). Particularly, the western Rif Mountain range is a real pole of diversification of Mediterranean species due to the high diversity of habitats ([12], [46]). This region represents the African side of the Mediterranean Intercontinental Biosphere Reserve between Morocco and Spain, which was designated by UNESCO’s Man and Biosphere Program in 2006. Except for tree species growing in arid or extremely cool conditions, such as Acacia Miller spp., Argania spinosa (L.) Skeels (argan tree) and Cupressus atlantica Gaussen (Atlas cypress) in southern Morocco or Juniperus thurifera L. (Spanish juniper) in the High Altas Mountains, all major Moroccan tree species are present in the western Rif. Some of them are endemic relict species such as Abies marocana Trabut (Moroccan fir) and Pinus nigra Arnold subsp. mauretanica (Maire & Peyerimh.) Heywood (black pine), while others are growing at the marginal areas of their distribution range, mainly Alnus glutinosa (L.) Gaertner (black alder) at its southern distribution limit and Cedrus atlantica (Endl.) Manetti ex Carriere (Atlas cedar) at its northwestern limit. Besides, many of these species are threatened, such as A. marocana and Cedrus atlantica, which are listed in the IUCN’s Red List as endangered ([4], [50]).

The Rifian forests, exposed to various historical influences, have developed a complex vegetation structures due to heterogeneous local ecological conditions and multiple forms of anthropic pressure ([49]). Indeed, in the western Rif, similar to many Mediterranean regions, the natural ecological balance, which was progressively established by the improvement of the climate since the beginning of the post-Würmian interglacial period, was profoundly affected by the rise of Neolithic civilizations that practiced pastoralism and agriculture ([46]). In this region, agriculture and livestock development began between 6,000 and 7,000 years ago, along with other practices such as clearing, fires, timber, and firewood collection. These activities intensified further in the 10^th century ([7]). During this period and given the rapid transformation, the effects of climate played only a secondary role in establishing and maintaining natural landscapes compared to the anthropic impact ([46], [7]). Except for specific mountain areas that have remained sheltered, it is challenging to get a precise idea of the original vegetation, especially in lowlands and deep soils. Nowadays, tree species (when still present) occupy only some parts of their original potential range. In this context, studies in the Rif on potential vegetation could only recognize and identify specific types based on natural vegetation fragments that are more or less preserved in relict habitats ([10]).

Faced with forest degradation stemming from multiple causes, commendable efforts have been made by the Moroccan State to conserve and restore forest stands over the past two decades, based on native species reforestation after having obtained, since the 1950s, very modest or even poor results using foreign species. In Morocco, designing forest restoration projects requires a thorough knowledge of local socio-economic aspects and a mastery of forestry techniques. The Moroccan government has launched an ambitious “Morocco’s Forests 2020-2030” strategy to address the socio-economic and environmental issues impeding forestry development. This strategy aims to manage the forestry sector sustainably and inclusively. A program of 600,000 ha of reforestation and rehabilitation in degraded ecosystems is part of the strategy’s action plan. The main focus of the strategy is to promote reforestation using species native to the country. However, choosing one or more target species that adapt to local ecological conditions is a prerequisite for any rehabilitation action and silvicultural intervention.

In this context, the current study utilizes species distribution models (SDM) as an effective method to identify and map the most suitable areas for a particular species within a specific habitat. We modeled the potential natural vegetation (PNV) at fine spatial resolutions using a predictive statistical GIS-based modeling approach (i.e., MaxEnt) to assist in decision-making regarding the selection of target species in management plans. MaxEnt is based on a maximum entropy algorithm ([44]) and was used to develop statistical predictive models for the potential distribution of fourteen Mediterranean trees native to Morocco. The modeling results are relative distribution probabilities for each species, and the species with the highest probability at a given site represents the potential species to be established at that site in reforestation activities.

SDMs, also known as ecological niche models or climate envelope models, are statistical predictive models based on environmental variables as explanatory variables to estimate species’ potential habitat ([28], [6]). BIOCLIM, the initial SDM package introduced in 1984, pioneered SDM techniques ([17]). SDMs are considered spatially continuous data on environmental suitability. They allow for the comparison of habitat factors present at species occurrence locations to those found within the study area, using various statistical methods ([35]). Ecological-niche characteristics based on the environmental conditions at known species presence sites allow a refined picture of their spatial distribution ([5]). For SDMs, the geographical coordinates of species occurrence represent a fundamental element of the input data, which consists of sampling locations found in scientific observations or collections under natural conditions, proving the species’ existence in a given region ([35]).

Key assumptions of SDMs are that species are at equilibrium with their environments and that relevant environmental gradients have been adequately sampled ([16], [24]). As a result, occurrence points theoretically reflect the habitat distribution that offers favorable conditions for a given species. In contrast, a lack of occurrence points may indicate that the habitat is unsuitable for the species or that the habitat has not been adequately sampled ([35]). The first reports of SDMs date back to the early 20^th century ([25]). The first documented correlative relationships between species and environmental variables date back to 1924 ([29]). Still, the statistical formalization of SDMs only began to take shape in the 1970s ([25]). Subsequently, in the 1980s, multivariate statistical techniques were developed with the BIOCLIM package ([17]). Since the 1990s, integrating GIS, statistical methods, machine learning, and increased data availability has further revolutionized SDMs ([25]).

Based on the maximum entropy theory ([44]), the MaxEnt model is recently becoming a common SDM used in ecology to identify the potential distribution of species due to its good predictive ability ([18], [45]). Compared with other models, MaxEnt appears preferable for its efficiency in making predictions over larger areas because it only requires presence data and fewer, easily obtainable environmental variables ([48]). Thus, it is considered among the best algorithms for providing accurate predictions ([24]).

In Morocco, applications of SDMs remain relatively scarce. Recently, MaxEnt modelling has been used in predicting the potential distribution of some species, such as A. marocana ([3]), Cedrus atlantica ([33]), and Quercus suber L. ([33]). While most of these studies focused on a single species, our study is the first that simultaneously examines fourteen forest tree species in Morocco, including seven broadleaved and seven coniferous trees, some of which are endemic. The broadleaved species consist of Quercus faginea Lam., Q. pyrenaica Willd., Q. rotundifolia Lam., Q. suber, Q. coccifera L., Q. lusitanica Lam., and Olea europaea var. sylvestris. Coniferous species are A. marocana, Cedrus atlantica, Pinus mauretanica, P. pinaster Aiton subsp. hamiltonii (Ten.) Huguet del Villar with two varieties (var. maghrebiana Huguet del Villar and var. iberica H. del Villar), P. halepensis Miller, and Tetraclinis articulata (Vahl) Masters. We modeled the potential distribution of these taxa under current ecological conditions and determined the key environmental variables that control their distributions. As these species are recognized as the main potential forest species in northwestern Morocco ([10]), a potential vegetation map model was also produced.

  Material and methods 

Study area

The present study was conducted in the Tangier-Tetouan region located in northwestern Morocco and known geographically as the western Rif (Fig. 1). It is located between latitudes 34° 30′ - 35° 55′ N and longitudes 06° 14′ - 04° 27′ W, covering an area of 12,650 km². It is bordered by the Atlantic Ocean to the west, the Gibraltar Strait to the north, and the Mediterranean Sea to the east. The highest mountain reaches up to 2159 m a.s.l. (i.e., Mt. Lakraa). The climate is Mediterranean with bioclimates ranging from semi-arid on the Mediterranean coast, subhumid to humid on the Atlantic coast and plains, to perhumid at high elevation. Based on the spatial distribution maps of bioclimatic parameters established by Ghallab & Taiqui ([27]), the mean annual precipitation is about 920 mm year^-1 with a maximum of 1830 mm year^-1 in high altitudes, and the mean annual temperature is about 17 °C, varying between 13 and 19 °C, with a minimum in January and a maximum in August.

Fig. 1 - Map of the study area and the location of the protected areas in northwestern Morocco.

  Enlarge/Shrink   Download   Full Width  Open in Viewer

According to the National Forest Inventory (NFI) statistics, the forest cover represents 32% of the study area ([30]). It is highly fragmented due to fires, clearings, overgrazing, and overexploitation. There are 19 Sites of Biological and Ecological Interest identified by the National Plan of Protected Areas carried out in 1996; they encompass different ecosystems and habitat types to be protected. These protected areas are highly rich in plant species which are utilized by humans in comparison to adjacent unprotected areas ([45]). Consequently, they are instrumental in biodiversity conservation and sustainable development.

Species occurrence data

The species studied and the number of occurrence sites used in our study are listed in Tab. 1. An overview of the biogeographical origin, bioclimates, altitudinal ranges ([9], [12]), and areas occupied (NFI stand type maps ([30]), revised based on recent forest management maps) of these species are also included in Tab. 1. All the studied species form metastable zonal forests in equilibrium with environmental conditions ([10]). Indeed, paleontological studies confirm the stability of the selected taxa in the study area. These taxa have been known since the basal Miocene and their importance increased in the Pliocene ([47]).

Georeferenced occurrence data has been collected from two main sources: (i) phytosociological studies by Benabid ([9], [10]) that mentioned sample locations; this list has been completed using GPS during our field visits; (ii) spatial database and maps from the NFI ([30]). In addition, we used the available maps from unpublished forest management studies and reports conducted by the Regional Directorate of the National Agency for Water and Forests Tangier-Tetouan Al Hoceima. Point extraction from these digital maps has been randomly generated using ArcGIS^® v. 9.3 (ESRI, Redlands, CA, USA) with a minimal distance of 50 m corresponding to the spatial resolution of our environmental database variables. Oversampling records in some areas were identified and removed by creating 50 m buffer zones around each observation point. A final dataset of 17,692 presence records, ranging from 74 to 9,004 occurrences per species, was included in the SDMs.

Environmental predictors

The explanatory variables are selected according to their strong predictive power of the potential distribution of the chosen species ([6]), their significant variability in the study area, and their independence or at least the lack of significant correlation with each other ([40]).

Based on investigations into the species’ ecology conducted in the study area ([10]), seventeen variables were initially identified as having the most significant influence on species habitat. The list included five climatic, one lithological, four topographic, four geographic variables, solar radiation, and two spectral indices (see Tab. S1 in Supplementary material).

The climatic variables were taken from the dataset compiled by Ghallab & Taiqui ([27]), who reported estimates from 1960 to 2000 with a spatial resolution of 50 m. This resolution was adopted for all remaining variables. The digital elevation model was obtained from the same database to derive the topographic variables and solar radiation. NDVI and TVDI were produced at a resolution of 30 m (i.e., the native resolution of the Landsat 8 image) and then reduced to 50 m using the nearest neighbor resampling technique in ArcGIS. Tab. S1 (Supplementary material) provides further details on how the selected variables were generated.

To reduce the multicollinearity among predictors ([44], [40]), we calculated the Pearson’s pairwise correlations between variables (Tab. S2), retaining only those that were not highly correlated (r < |0.70| - [20]). In the case of high correlation, factors with indirect ecological influence were removed ([28]), such as topographic and geographic variables in our study. Since TVDI (i.e., a measure of soil moisture status) is mainly developed based on NDVI, which is directly related to vegetation types, NDVI was eliminated even though the correlation between the two indices was relatively low. Overall, ten variables were retained for the analyses: mean annual precipitation, mean annual temperature, mean minimal temperature of the coolest month, mean maximal temperature of the hottest month, slope, curvature, aspect, lithology, solar radiation, soil humidity, and TVDI.

Modeling approach

Due to the lack of data on the biology of some studied species and considering the spatial dimension of the study area, the correlative approach proved to be more appropriate than the mechanistic approach for determining the ecological niche of species ([28]). Based on the comparative examination of modeling methods ([24]), we used the MaxEnt v. 3.3.3 software to model the current distribution of selected species. This method makes it possible to predict a probability of presence between 0 and 1 from floristic records and current ecological variables ([28]). The prediction maps can be imported into a GIS software for further analysis. The maps were created at the resolution of the original predictor datasets (50 m).

The modeling was run for each species to build a complete model encompassing all ten predictor variables. This enabled the determination and comparison of the importance of each variable’s contribution to the distribution of each target species. For this purpose, MaxEnt’s jackknife test was used ([44]).

The generated models were calibrated employing 75% of the occurrence records for each species as training data, while the remaining 25% served as test data for model validation. The Bootstrap method was applied with 10 iterations. A default regularization multiplier of 1 was applied to reduce overfitting, and a maximum iteration number of 1000 was chosen to ensure model convergence. All other parameters have been set to the recommended default values.

The performance of the model to predict the potential distribution was assessed by calculating the AUC (Area Under Curve) with the ROC (Receiver Operating Characteristic) curve ([44]). From the ROC curve, the AUC value is calculated to test the relevance of the model. AUC takes values between 0.5 and 1. According to Hosmer & Lemeshow ([31]), a value greater than 0.9 corresponds to an exceptional model, indicating that the adjustment of the model is perfect; a value between 0.9 and 0.8 assumes an excellent model, a value between 0.8 and 0.7 characterizes an acceptable simulation, and a random model is defined by a value of 0.5 which means that the model has no discriminant value. The obtained results of the current SDMs were exported into ArcGIS.

To perform statistics on the predicted SDMs compared to the observed ones (NFI map of 2014), binary maps were developed by segmenting the predicted areas of each species into two classes (suitable/unsuitable) based on threshold values of maximum training sensitivity plus specificity ([36]). Thus, for a probability of occurrence lower than the value of this threshold, the habitat was qualified as unfavorable; otherwise, the habitat was considered favorable.

We used the “HighestPosition” tool (Spatial Analyst) in ArcGIS to generate the PNV map. This tool allows to determine the species with the highest occurrence value for each cell by analyzing a set of species occurrence rasters. As a result, it identifies the raster corresponding to the species that is prevalent in each specific location. To assess the agreement between the predicted PNV map and the current species distribution, we used the Kappa (K) statistic, which is widely used to test the reliability of vegetation maps. The K value generally ranges from 0 to 1, and to categorize the agreement, the following classes are used: very poor to poor if K < 0.4, fair if 0.4 < K < 0.55, good if 0.55 < K < 0.7, very good if 0.7 < K < 0.85, and excellent K > 0.85 ([42]). This test was performed using SPSS^® v. 20 (IBM, Chicago, IL, USA) based on the information extracted from a sample of 5605 points systematically distributed over the study area. The workflow of our modeling is summarized in Fig. 2.

Fig. 2 - Workflow diagram.

  Enlarge/Shrink   Download   Full Width  Open in Viewer

  Results 

Model accuracy and contribution of environmental variables

The quality of the SDMs developed for the fourteen species was considered exceptional for ten species, excellent for one species, and acceptable for the remaining three, as none of the models had an AUC less than 0.70 (Fig. S1, Fig. S2 - Supplementary material).

The relative contributions of the environmental variables to the construction of SDMs are presented in Tab. S3 (Supplementary material). Abies marocana distribution is mainly structured by the variables in the following order: Tmoy, Tmax, Pmoy, soil humidity and lithology. Tmoy and Pmoy contributed with the highest percentages to C. atlantica spatial distribution followed by Tmax, lithology, and TVDI. As for A. marocana, the Tmoy is the variable significantly affecting cedar distribution model. Similarly, the Tmoy is the most explanatory variable for the distribution of P. maghrebiana, contributing 81.1%. Lithology, followed by Tmax, TVDI, and Pmoy, fills in the rest of the model construction.

The distribution of O. sylvestris is primarily and equally driven by Tmoy and Pmoy, followed by Tmax. The distribution of P. halepensis is essentially determined by the following variables (in descending order of importance): Tmoy (48%), slope (17%), Tmax (13%), lithology, exposure, and soil humidity. In the case of P. iberica, Tmax was the variable that contributed the most (63.6%), followed by Tmin, TVDI, slope, Tmoy, lithology, and Pmoy. For P. mauretanica, the Tmoy was the only variable explaining its distribution (97.7%) and the most contributing to the model construction (75.6%). Lithology, Tmax, TVDI, and Tmin were the other significant predictors (Tab. S3 in Supplementary material).

Soil humidity contributed by 52.7% to the explanation of Q. faginea distribution, while climatic variables (Tmax, Tmoy, and Pmoy) contribute to the rest. Soil moisture also had the highest weight in the construction of the model. The distribution of Q. coccifera was primarily related to climatic variables (Tmax, Tmoy, and Pmoy) with a total weight of 75%, followed by slope, solar radiation, Tmin, and TVDI. These variables retain the same order of importance in explaining the species distribution or constructing the model. The main variables in driving the Q. lusitanica distribution were Tmax, Tmoy, lithology, Pmoy, soil humidity, and Tmin. For Q. pyrenaica, the main variables were Tmoy, Tmin, TVDI, Tmax, Pmoy, and slope (Tab. S3), while the Tmin had the most significant weight in the construction of the model. For Q. rotundifolia, Tmoy was the main driving factor of the distribution (71.6 %) and construction of the model (81.3%), followed by Tmax and slope. The distribution of Q. suber was largely determined by lithology (63.5%) followed by Pmoy, then Tmax and TVDI; these variables are also arranged according to this same order of importance to build the model. Finally, for T. articulata, climatic variables (Tmax, Tmoy, and Pmoy) and slope were the main drivers of the distribution (96%) and the model construction (86%).

The results of the jackknife test (see Fig. S3 and Appendix 1 in Supplementary material) showed that excluding individual variables did not lead to a significant decrease in training gain, with the exception for several species such as P. mauretanica, Q. lusitanica, and Q. suber. This finding implies that no single variable provides additional valuable information beyond what was provided by the other variables. Moreover, for all species analyzed, none of the models developed by omitting any variable and using the remaining variables exceeded the performance of the models that incorporated the whole set of variables. This finding underscores the critical importance of including the entire array of variables to enhance model efficiency in the study area.

Mapping of the species potential ranges and comparison with current distribution

The predicted distribution of A. marocana is much larger than its current distribution observed (Fig. 3, see also Tab. S4 in Supplementary material). The predicted potential area extends over 22,429 ha and over 18,920 ha, excluding the unforested lands. The mapped area is 3,000 ha, i.e., 13% and 16% of the predicted potential area with and without unforested lands, respectively. C. atlantica potential area is discontinuous and distributed over the main mountain summits as patches of unequal extension. The predicted potential area is 36,198 ha, of which natural spaces still occupy 29.736 ha. The mapped area is 4,192 ha, i.e., 12% and 14% of the total predicted potential area with and without unforested area, respectively (Fig. 3, see also Tab. S4). Pinus maghrebiana clearly shows high ecological plasticity. While its current range covers only 9,597 ha, the predicted potential area, including unforested land, is 60,727 ha. The species’ current distribution represents 16% and 20% of the potential area with and without unforested lands, respectively (Fig. 3, see also Tab. S4 in Supplementary material)

Fig. 3 - Superposition of the current distribution (red color) and the predicted potential area (green color) of the fourteen species studied.

  Enlarge/Shrink   Download   Full Width  Open in Viewer

The suitable area predicted for O. sylvestris is large in the study area, equalling 326,562 ha and only 57,876 ha if unforested lands are not considered. While the current range of P. halepensis covers only 1,382 ha, the predicted potential area extends over 65,834 ha including unforested lands. In the case of P. iberica, the potential area was discontinuous and theoretically extends over 66,285 ha. The actual distribution covers only 1,229 ha, i.e., 2% and 3% of the potential area without and with unforested area, respectively. P. mauretanica potential area extends over 10,658 ha, with 9,623 ha outside the unforested areas, while it currently covers only 501 ha, i.e., about 5% of its potential area.

Quercus faginea showed a total predicted potential area of about 54,279 ha, while the mapped range was about 10,790 ha, i.e., 20% of the potential area. Q. coccifera had the largest potential area among all studied species, extending over 351,736 ha, including unforested lands, and only 127.020 ha if excluded. The degradation of Q. coccifera stands seems clear, as the existing map only identifies 17,722 ha, equalling 5% of the total predicted area (Fig. 3, see also Tab. S4 in Supplementary material) or 14% if unforested areas are excluded. The total potential area of the Q. lusitanica extends over 65,254 ha and only 35,847 ha if the unforested lands are excluded. It should be noted that this is the first mapping carried out for Q. lusitanica. The existing cartography confirmed that the predicted distribution of Q. pyrenaica corresponds to its known range. Indeed, the predicted potential area was 45.744 ha, while the area occupied is only 2,087 ha (5% of its total potential area). The total predicted area of Q. rotundifolia was large and extended over of 238.457 ha, and only 124.481 ha if bare soil is excluded. The stands currently mapped concern only 42,336 ha, i.e., 18% of the total potential area and 34% of the silvicultural areas. The predicted favorable area for Q. suber covers an area of 316,874 ha against 173,641 ha if bare soil is excluded. The occupied area was about 113,330 ha, i.e., 36% of the potential area and 65% of the zones currently covered by natural woody stands.

The suitable area for T. articulata extended over a total potential area of 188,800 ha, while it is reduced to only 76,071 ha if bare land is excluded, which indicates the high degradation this ecosystem has experienced. The area currently occupied by T. articulata is about 21,205 ha, being 11% of the total predicted potential area and 28% if bare land is excluded (Fig. 3, Tab. S4 in Supplementary material).

The loss of potential forest areas due to their conversion into non-forested lands had varied effects on the species studied. Notably, O. sylvestris is the most impacted by deforestation, with 82% of its potential area converted to agricultural or urban use. Q. coccifera ranked second with a 64% loss of area, while T. articulata experienced a 60% loss. The remaining species had suffered the following losses: Q. rotundifolia (-48%), Q. suber (-45%), Q. lusitanica (-45%), P. iberica (-36%), P. halepensis (-22%), P. maghrebiana (-21%), C. atlantica (-18%), A. marocana (-16%), Q. pyrenaica (-14%), Q. faginea (-11%), and P. Mauretanica (-10%).

Species distribution in the environmental

The 2-dimensional environmental space defined by the Pmoy and Tmoy variables shows that most of the study area falls in the environmental gradients between 600 mm to 1200 mm and 16.5 °C and 19 °C (Fig. 4). The relief, especially the high altitudes, and the distance from the sea coast increase the size of this space. Most broadleaved trees tend to occupy places with moderate Tmoy and high Pmoy (Fig. 4a). Conversely, conifers prefer relatively cold sites with sufficiently abundant water balance (Fig. 4b).

Fig. 4 - Potential areas of species (a: broadleaved; b: conifers) in the two-dimensional environmental space formed by mean annual precipitations (Pmoy) and mean annual temperatures (Tmoy). Species labels - Qp: Q. pyrenaica, Qc: Q. faginea, Qr: Q. rotundifolia, Qs: Q. suber, Ql: Q. lusitanica, Qcc: Q. coccifera, Oe: O. sylvestris, Am: A. marocana, Ca: C. atlantica, Pn: P. mauretanica, Pmag: P. maghrebiana, Ph: P. halepensis, Pib: P. iberica, Ta: T. articulata.

  Enlarge/Shrink   Download   Full Width  Open in Viewer

The relict species, namely A. marocana, C. atlantica, and P. mauretanica, are confined in a limited environmental gradient, i.e., the highest altitudinal zone where bioclimate is the wettest and coldest. Tmoy ranges for the three species are 1.7 °C, 2.6 °C, and 1.5 °C, respectively. This may indicate that these species would be extremely sensitive to global warming. The ranges of the ten environmental variables that correspond to high (≥50%) occurrences for each taxon are shown in Tab. 2.

Mapping of the potential natural vegetation (PNV)

The map of PNV is represented in Fig. 5a. The results of the binomial test of the comparison between this map obtained by the species prevalence domain model and the current mapped species distribution (NFI map of 2014) are reported in Tab. S5 (Supplementary material). The test showed that out of 5,605 sampled points, 4,282 (76%) agreed between the predicted and mapped species types. Using a theoretical probability (proportion test) set at 76%, the p-value obtained (0.249) is greater than 0.05 (significance level). On the other hand, the Kappa statistic obtained by comparing the two maps gives an overall value of 0.67, indicating no significant difference between the two maps. The visual comparison of our simulated PNV map (Fig. 5a) with that developed by Benabid ([10]) based on phytosociological studies (Fig. 5b) revealed that the former provides valuable detail on the distribution of lowland species that have historically undergone the most extensive clearing, thereby constraining their distribution. Compared to the map of Benabid ([10]), Q. coccifera would potentially occupy marl soils over a much larger area. In addition, O. sylvestris could potentially extend its range to clay soils in the plains, valleys and lowlands in the study area. T. articulata and sclerophyllous broadleaf species (i.e., O. sylvestris, Q. coccifera, Q. suber, Q. rotundifolia) are the only species whose potential range in the PNV map exceed 50% of their individually modelled potential areas (Fig. 5).

Fig. 5 - Potential Natural Vegetation (PNV) map of the western Moroccan Rif simulated under current environmental conditions.

  Enlarge/Shrink   Download   Full Width  Open in Viewer

  Discussion 

Based on SDMs, we predicted for the first time the potential distribution area of fourteen major Mediterranean tree species in the western Rif of Morocco, emphasizing the main environmental factors that drive their potential distributions. These species were selected based on the historical value in the region and their stability under the local ecological conditions. The dataset of the investigated environmental variables included ten factors that may influence the distribution of tree species in the study area (Tab. S2, Tab. S1 - Supplementary material). Based on the AUC values, the predicted distribution maps of different species were very satisfactory. The most relevant models concern endemic species or specialist species with limited distribution areas due to their narrow ecological niche ([38]), such as A. marocana, P. mauretanica, P. maghrebiana, C. atlantica, Q. pyrenaica, Q. lusitanica, and P. iberica.

The SDMs allowed modeling the ecological behavior a given species along the gradients of the ten explanatory variables. Overall, the bioclimatic factors (Tmoy, Tmax, Pmoy, and Tmin) had the highest contribution in SDMs, followed by soil variables (lithology and TVDI), topographic variables (slope, aspect and curvature), and finally luminosity (solar radiation), with 76%, 17%, 7%, 1%, respectively (Tab. S3 in Supplementary material). Climatic variables (mainly precipitations and temperatures) represent the most driving factors of SDMs, since each species can tolerate a certain threshold of these factors ([37]).

Tmoy proved to be the primary discriminant factor for many of the species analyzed, while Pmoy generally played a secondary role in SDMs, despite it significantly affects several species such as O. sylvestris, C. atlantica, A. marocana, and P. maghrebiana. Species growing in the Tingitane peninsula or the Mediterranean side (P. iberica, Q. coccifera, and Q. lusitanica) were greatly influenced by Tmax, which is significantly mitigated by Mediterranean Sea effects. The ecology of Q. faginea and Q. suber is linked to soil characteristics; the first depends on more humid soil while the second is mainly conditioned by the nature of the substrate. Q. pyrenaica distribution seems to be influenced by Tmin and the contributions of soil humidity, though Pmoy and Tmoy were also important for this species. Indeed, it seems that the ecological niches of the studied species are mainly linked to bioclimatic factors, especially temperature. These results agree with previous findings by Bravo et al. ([19]), who stated that temperature is the predominant ecological factor in the Mediterranean mountains. However, water also plays an important role, making these ecosystems vulnerable to climate warming. Based on our results, the defined ranges of the ten environmental variables for each species showed significant consistency with the ecological requirements of these species.

The following section discusses the findings for each of the studied species.

A. marocana: The species grows on the wettest slopes of the Moroccan most humid region, on a calcareous-dolomitic substrate between 1400 and 2100 m of altitude with a mean annual rainfall ranging from 1940 to around 2000 mm ([9]). Our results support the previous knowledge of this endemic species and further clarified its bioclimatic environment. Indeed, the presence probability of A. marocana is the highest where Tmoy is lower than 13.5 °C, Tmax ranges 31.5-32.5 °C, and Pmoy is higher than 1560 mm (Tab. 2, Fig. 3). The potential area of A. marocana was 22.429 ha including unforested areas but 18.920 ha if these are excluded. These extensions are wider than those obtaines by Alaoui et al. ([3]), who found only 13.084 ha. For other Mediterranean firs, López-Tirado & Hidalgo ([37]) found that A. pinsapo is dependent on the winter rainfall and the mean coldest temperature.

C. atlantica: The most suitable habitat for this species is observed in the Rif between 1400 and 2300 m of altitude on various substrates ([9]). However, siliceous bedrocks (i.e., schist and sandstones of the Rif, basalt of the central Middle Atlas) seem to offer a more favorable water balance than dolomitic limestone substrates (eastern Middle Atlas, eastern High Atlas - [11]). The Pmoy requirements of C. atlantica are between 1390 and 1786 mm in the western Rif ([41]), consistent with our results. The presence probability of C. atlantica is higher where Tmoy is less than 13.7 °C, Tmax ranges 30.5-32.5 °C, and Pmoy is higher than 1400 mm (Tab. 2, Fig. 3).

P. iberica: It is an Iberian-Rifian endemic variety of P. pinaster. In Morocco, it only exists in the eastern part of the Tingitane peninsula ([10]). Given the limited distribution of this variety, there is a lack of knowledge of its ecology. Moreover, no action was taken to improve P. iberica stands, although it is recommended for reforestation programs at the regional scale. Based on our results, P. iberica can extend over a much larger area than its current range (Fig. 4). It showed high plasticity in regards to the temperature and precipitation, where Tmin, Tmax, and Pmoy are 5.5-8.5 °C, 29.5-32.5 °C, and 180-1200 mm, respectively (Tab. 2, Fig. 4). These results support the potentialities of P. iberica as an important plantation species in northern Morocco ([13]).

P. maghrebiana: It is a Moroccan endemic variety of P. pinaster. Our results represent a complement to its existing ecological knowledge. P. maghrebiana occupies the E, NE, or SE slopes at an altitudinal range of 700-1700 m, and it is indifferent to the substrate nature, although more present on limestone substrates ([10]). The high probabilities of P. maghrebiana occurrence are found where Tmoy is 12-14 °C, Tmax is 31-32 °C, on carbonate or siliceous substrates (Tab. 2). Like P. iberica, P. maghrebiana showed high plasticity in regards to the temperature (Fig. 4), but it is more constrained by precipitation (> 1200 mm), which may limit its role in reforestation to some humid sites of northern Morocco ([13]).

P. mauretanica: This species grows in very small clumps only on E, NE, and SE slopes of the eastern part of the Rif calcareous ridge summits, between 1300 and 1700 m a.s.l. ([10]). Notably, the most discriminant variable correlated to the high probability of its presence was Tmoy (12.7-13.3 °C), while other variables seem unrelated, though Pmoy was high (1450-1800 mm). However, in southern Spain, winter rainfall and mean coldest temperature have been reported to be the most influencing factors on the distribution of Pinus nigra subsp. salzmannii ([37]). The mapping of P. mauretanica clearly showed that the high probabilities of occurrence are limited to a very small area, suggesting that it has a specialist ecological behavior and that the study area only offers a range close to favorable environmental conditions. Our results represent a useful complement to the ecological knowledge of this species, which was limited so far.

P. halepensis: It forms scattered stands at plains, low, and mid-mountains from the semi-arid to subhumid bioclimates on all soil types (though rarely on siliceous soils), with the best stands located on limestone marl ([15], [10], [26]). P. halepensis is distributed along the Mediterranean coast as stunted clusters under the sea spray effect. In the western Rif, it is located on ridges or E and SE slopes ([10]). While Pmoy and Tmin appear to have insignificant influence, we found higher probabilities of occurrence under Tmoy of 15.5-18 °C and Tmax between 32.2 and 34.5 °C; these favorable climatic conditions are mainly localized in the east of the study area. These results support the available knowledge on the species, as its distribution is reported to be driven by relatively higher mean temperature and lower rainfall (around 250-450 mm - [15]). Interestingly, on the NE coast of Jbel El Kebir (western Tangier), which was predicted to be suitable for the species, we located many individuals, most of which had a lying shape due to the almost permanent winds. The reduction of P. halepensis area may be due to its low resistance to competition with Q. rotundifolia which strongly resprouts after fires or cuttings ([8]). Indeed, phytosociological studies considered that most P. halepensis stands represent transient or paraclimax groups whose normal evolution leads to other climax forests, particularly Q. rotundifolia ([2], [10]).

T. articulata: It is an endemic North African conifer, but it is in Morocco that it grows better, forming several associations ([2], [9]). As a thermophilic and xerophilic species, T. articulata occupies an altitudinal range between 0 and 500 m in the western Rif, and it is essentially confined to limestone and calcareous-marly substrates ([9]). In the study area, we found that T. articulata had a higher probability of occurrence where Tmoy range is 14.4-16.5 °C, Tmax 29.2-32.4 °C, Tmin 6-11 °C, and Pmoy < 800 mm. These results align with the available knowledge on the species, i.e., Pmoy of 600-950 mm, and Tmin of 5.5-8°C ([10]). However, the species spans over larger ecological gradients than those reported above, as its range extends outside the study area.

O. sylvestris: It grows in the study area between 0 and 500 m a.s.l., with precipitation around 600-950 mm ([9]). It is indifferent to soil nature and is one of the few species that can grow on clay soils, although it prefers limestone and schist ([14]). Currently, the species climax has almost disappeared, and the remaining populations cover rocky steep soils unsuitable for cereal crops or arboriculture ([10]). The occurrence of O. sylvestris has higher probability where Tmoy, Tmax, and Pmoy are in the following ranges: 17.5-20°C, 31.5-37.2°C, and 700-950 mm, respectively (Tab. 2), which are consistent with the available ecological description of the species. The potential predicted area of O. sylvestris was the largest among the fourteen studied species (Fig. 4). However, 82% of this area has been converted into croplands or arboriculture (mainly olive trees). Thus, this suggests that O. sylvestris underwent the most extensive deforestation throughout the history of human settlement in the region.

Q. faginea: In Morocco, this species is represented by three subspecies, all present in the Rif ([10], [26]): ssp. canariensis Willd. grows on siliceous substrate at an altitude range of 130-800 m in the Tingitane peninsula and 900-1700 m in the western and central Rif; subsp. faginea grows on limestone and in the west of Rif at 1400-1700 m of altitude; and subsp. telemcensis Mair et Weiller in the Rif. Generally, Q. faginea is indifferent to the soil nature, and its diffusion is mainly related to bioclimatic conditions ([2]). However, the species colonizes deep, slightly sloping colluvial soils at the foot of cliffs and cornices or in cool valleys ([10]) generally charaterized by a good water balance. We found that the high probabilities of Q. faginea occurrence were determined by temperature and precipitation (Tab. 2), which perfectly matched the available knowledge on the species ecology. Its stands appear relatively conserved, with only 11% of the potential area being converted to croplands (Tab. S4 in Supplementary materials).

Q. coccifera: It is found in the Rif, Beni-Snassen, and Taza, between 0 and 1000 m in the western Rif on limestone substrates, marl-limestone, rarely schist or sandstone, under perhumid and humid bioclimates ([9]). The species’ climax is only observed in marabouts ([10]). The available mapping of the species only delimits high-altitude stands while neglecting the lower ones. Generally, Q. coccifera pure stands are rarely observed and often considered matorrals. We found high probabilities of its occurrence under Tmoy of about 17.7-19 °C, Tmax 28.8-36.5 °C, Tmin 4-11 °C, and Pmoy 800-1100 mm (Tab. 2).

Q. lusitanica: It is an Iberian-Rifian endemic species located in the western and southwestern Iberia and the Tangitane peninsula ([26]). It grows between 100 and 500 m of altitude, exclusively on sandstone, and forms low matorrals (< 1.5 m of height - [10]). It benefits from maritime influences, which ensure intense atmospheric humidity, promoting a reduction in temperature variations. It should be noted that no mapping exists for Q. lusitanica, thus, we provide the first species distribution map. Its occurrence seems to be relatively influenced by similar temperature ranges (Tmoy and Tmax) compared to Q. coccifera, while Tmin and Pmoy seem to have insignificant effects (Tab. 2).

Q. pyrenaica: It thrives on the siliceous western upper Rif, characterized by heavy precipitation, growing on moderately deep but regularly cool sloping soils at 1200-1700 m of altitude ([10]). In contrast to the previous Quercus species, Q. pyrenaica requires a specific habitat, with a high probability of occurrence where Tmoy is lower than 15.5 °C, Tmin between -1 and -2 °C, and Tmax of 33.9-34.5 °C (Tab. 2).

Q. rotundifolia: Covering an area of about 1.342.000 ha, it is the most important forest species in Morocco ([2]). It appears to be the most ecologically plastic species which grows under many varied and harsh conditions. It can withstand a variation in Tmin ranging from -3 °C to +7 °C ([14]). However, we found that the distribution of this oak was insensitive to Tmin and that is limited to a Tmax not exceeding 35.5 °C, in contrast to a value of 42 °C reported by Dahmani ([23]). Q. rotundifolia seems to grow on all types of substrates, except on compact, anoxic, or saturated soils ([1]).

Q. suber: It is an endemic of the Atlantic side of the Mediterranean basin ([47]). In Morocco, Q. suber forests extend from plains to mid-mountains (0 to 1600 m) in the Rif, eastern Middle Atlas, High Atlas, and in the Atlantic plains from Rabat to Casablanca ([10]). It is a light-demanding tree growing exclusively on siliceous soils (i.e., calcifuge - [2]). In contrast to Benabid ([10]), who reported that Q. suber requires rainfall ranging from 650-700 to 1700 mm, and Tmin from -2°C to 10°C, we found the highest probability of occurrence where Pmoy is 1000-1650 mm and Tmin 0.5-5 °C (Tab. 2).

Regarding the six Quercus species, we found that temperature and precipitation are the most critical environmental variables driving their occurrence. Generally, they require a high amount of precipitation, e.g., more than 1400 mm for Q. faginea and Q. pyrenaica, and more than 1200 for Q. Rotundifolia (Tab. 2). Based on Fig. 4, Q. suber and Q. rotundifolia appear to be the most favored species; Q. suber supports warmer conditions with no less than 400 mm of Pmoy, while Q. rotundifolia better tolerates cold conditions and precipitation variations. This last feature is partially shared with P. halepensis. Current conditions allow the relict A. marocana, C. atlantica, and P. mauretanica to dominate only in restricted areas.

In the present study, we modeled the northwestern Morocco PNV formed by the prevalence domains of major woody species. The overall accuracy (K = 67%) was similar or slightly better than those obtained in previous studies (e.g., 64% in [34], 58% in [43], and 48% in [32]).

Conservation and restoration implication

The maps produced delineate the potentially suitable areas for fourteen Mediterranean tree species in the study area. The calibration of the models is based on the occurrence sites in the current species distribution area, which is determined by complex interactions of biotic and abiotic factors. Based on the ecological niche principle, spatial distribution modeling improved our understanding of the ecological behavior and the potential range of species, particularly those endangered such as A. marocana, C. atlantica, and P. mauretanica. This undoubtedly represents a fundamental tool for implementing adequate conservation strategies. Moreover, our maps created by the SDM approach, whether for individual species distribution or integrated PNV maps, would strongly contribute to prioritizing intervention sites and conceiving restoration and rehabilitation programs for degraded areas using the appropriate species. The species studied currently cover only small parts of their potential ranges due to historical deforestation and interspecific interactions. The PNV model maps the species distribution with satisfactory accuracy (76%) due to these interactions (realized niche).

Nowadays, ecosystem restoration is a top priority worldwide and is on the agendas of many international, regional, and national organizations, such as the World-Wide Fund for Nature (WWF), the Food and Agriculture Organization of the United Nations (FAO), the European Commission, the Society for Ecological Restoration (SER), and the United Nations Environment Program (UNEP). The alarming status of ecosystem degradation threatens Africa’s ability to meet the sustainable development goals outlined in the 2030 Agenda, the African Union’s Agenda 2063, and other initiatives ([21]). To stop forest degradation and restore degraded forests, the African Union has created a comprehensive guideline, “Sustainable Forest Management Framework for Africa (2020-2030)”, which may be used as a reference for these countries and other African stakeholders involved in the forestry sector.

Our findings should be complemented by modeling the potential distribution of the studied species under various climate change scenarios. This will help in formulating urgent strategies for the sustainable management of forest ecosystems. Additionally, future studies should compare multiple modeling techniques to reduce biases and more accurately predict the current and future distribution of species. In this context, the mechanistic approach, based on theoretical or experimental knowledge of the relationship between species physiology and abiotic factors, has gained significant interest among modelers.

  Conclusions 

In the present study, SDMs allowed the prediction of the potential distribution areas of fourteen Mediterranean tree species in the western Rif (Morocco) with a high degree of accuracy. MaxEnt revealed high predictive power for ten species, mainly specialist and endemic, and an excellent to acceptable power for the four remaining species. Our results provide detailed information on species ecology by examining their ecological behavior along the gradients of ten variables, thus better defining the ecological niche of threatened species such as A. marocana, C. atlantica, P. Mauretanica., and (for the first time) Q. lusitanica. Moreover, our findings suggest the primary role of climatic factors (Tmoy, Tmax, Pmoy, and Tmin) in the delimitation of species’ ecological niches. The potential distribution area for all fourteen species was larger than the current area. Species in plains or mid-mountains were the most affected by deforestation, as a consequence of agricultural expansion. For example, 82% of the O. sylvestris potential area was affected by past clearings, followed by Q. coccifera, and C. atlantica, which lost 64% and 60% of their areas, respectively.

The relict species A. marocana, C. atlantica, and P. mauretanica are confined to the highest altitudes, making them highly vulnerable to global warming. Therefore, the long-term temperature goal of the Paris Agreement (a 1.5 °C rise in global temperatures) would not be sufficient to ensure the survival of these species.

  List of abbreviations 

MaxEnt: Maximum Entropy, PNV: potential natural vegetation, SDMs: Species Distribution Models, GIS: Geographic Information System, NFI: National Forest Inventory, NDVI: Normalized Difference Vegetation Index, Pmoy: mean annual precipitation, Tmax: maximum temperatures of the hottest month, Tmin: minimum temperatures of the coldest month, Tmoy: mean annual temperatures, TVDI: Temperature Vegetation Dryness Index, IUCN: International Union for Conservation of Nature, K: Kappa statistic.

  Acknowledgments 

AG conceived the study and the methodology, conducted the data analysis, wrote the original draft, and supervised the manuscript. HB and SL contributed to the conceptualization, the result analysis, and the writing and revising of the manuscript. MB helped to draft and revise the manuscript.

We thank the Regional Directorate of the National Agency of Water and Forests Tangier-Tetouan Al Hoceima for providing the maps of the National Forest Inventory (NFI). We express our gratitude to Prof. Lahcen Taïqui for his valuable insights and assistance in conducting this research. We are grateful to the four anonymous reviewers for their helpful comments that improved the quality of the manuscript.

Funding: this research did not receive any funding.

  References

  Supplementary Material

Tab. S1 - The seventeen environmental variables used in the species distribution modeling in the present study.
Tab. S2 - Pearson correlation matrix for the ten variables selected in the modeling process.
Tab. S3 - Relative contributions (%) of environmental variables in constructing species distribution models.
Tab. S4 - Summary of the predicted potential and current areas of the fourteen taxa studied. The current area is based on the National Forest Inventory (NFI).
Tab. S5 - Results of the binomial test and the Kappa statistical index obtained after comparing the developed PNV model with the real distribution of the studied species.
Fig. S1 - Omission rate and predicted area for broadleaved and coniferous species.
Fig. S2 - Receiver operating characteristic (ROC) curve for broadleaved and coniferous species.
Fig. S3 -Results of the jackknife test showing the importance of the ten studied variables in determining the species potential range using training gain.
Appendix 1 - Rationale of the jackknife test applied.