The Relict Ecosystem of Maytenus senegalensis subsp. europaea in an Agricultural Landscape: Past, Present and Future Scenarios

: Maytenus senegalensis subsp. europaea is a shrub belonging to the Celastraceae family, whose only European populations are distributed discontinuously along the south-eastern coast of the Iberian Peninsula, forming plant communities with great ecological value, unique in Europe. As it is an endangered species that makes up plant communities with great palaeoecological signiﬁcance, the development of species distribution models is of major interest under different climatic scenarios, past, present and future, based on the fact that the climate could play a relevant role in the distribution of this species, as well as in the conformation of the communities in which it is integrated. Palaeoecological models were generated for the Maximum Interglacial, Last Maximum Glacial and Middle Holocene periods. The results obtained showed that the widest distribution of this species, and the maximum suitability of its habitat, occurred during the Last Glacial Maximum, when the temperatures of the peninsular southeast were not as contrasting as those of the rest of the European continent and were favored by higher rainfall. Under these conditions, large territories could act as shelters during the glacial period, a hypothesis reﬂected in the model’s results for this period, which exhibit a further expansion of M. europaea ’ s ecological niche. The future projection of models in around 2070, for four Representative Concentration Pathways according to the ﬁfth report of the Intergovernmental Panel on Climate Change, showed that the most favorable areas for this species would be Campo de Dal í as (southern portion of Almer í a province) as it presents the bioclimatic characteristics of greater adjustment to M. europaea ’ s ecological niche model. Currently, some of the largest specimens of the species survive in the agricultural landscapes in the southern Spain. These areas are almost totally destroyed and heavily altered by intensive agriculture greenhouses, also causing a severe fragmentation of the habitat, which implies a prospective extinction scenario in the near future. concentration pathways (RCP2.6, RCP4.5, RCP6.0, RCP8.5); compare the retrospective and prospective results with the ecological niche model obtained from bioclimatic current data.


Study Area
The study area is located in the south and southeast of Spain. More specifically, it corresponds to the coastal area biogeographically encompassed by the Baetic and Murcian-Almeriensian provinces. This territory is characterized by a Mediterranean climate, with thermo-Mediterranean thermotype, modulated from dry to semiarid ombrotype, according to the bioclimatic classification proposed by Rivas-Martínez et al. [54]. It coincides with hotspots for plant biodiversity in southern Spain in terms of rarity and endemicity [55][56][57], and includes semiarid, as well as arid areas, which are considered among the most vulnerable ecosystems to global change drivers [58][59][60].

Species Distribution Models
Species distribution models (hereinafter SDMs) have been globally recognized as a useful tool in nature conservation and management, for instance, to refine the threat status of a species [61][62][63][64]. When applied to distribution data, they can predict distributions across geographic landscapes by multiple responses, improve image analysis or remote-sensing in order to lead the search for poorly known species [65][66][67][68], thus providing perceptions into the species' habitat, range and abundance [69][70][71][72][73]. Furthermore, several authors like Elith and Leathwick [74], Benito et al. [75], Fois et al. [76], and De Luis et al. [77,78] used SDMs based on the extant localities, as well as the respective current and future climate scenarios to predict the possible variation in the environmental niche of certain plant species, inferring ecological and evolutionary insights. MaxEnt 3.4.1 [79] was used to model the potential habitat of M. europaea in the different proposed scenarios. The application of the Maximum Entropy principle to estimate ecological niche modelling and potential distribution areas follows the studies of Phillips et al. [80] and Phillips and Dudík [81]. Said software has been used in many fields of science which proved its validity [81][82][83][84][85], being widely recognized as the most used tool, even for small sample sizes and poorly known species' distributions [71]. The modelling process is able to indicate suitable conditions for species in areas where their presence has not been registered or evidenced; several studies supported the reliability of SDMs [67,86,87]. Since this process is iterative, the modelling outcomes can highlight other factors related to the physical environment, anthropogenic influences or conditions of growth and reproduction of the species [88][89][90].
MaxEnt can operate with information about presence, which often represents the largest set of available data [91]. Only-presence-data algorithms usually represent the spatial distribution of the fundamental ecological niche of a species, while algorithms based on presence-absence data characterize more closely the distribution of the effective ecological niche [92]. The present study has been conducted with presence data alone, since a presence-absence modelling might generate controversies in light of the recent dynamics of the largest M. europaea population. Land changes due to the evolution of intensive crops in the province of Almería [25,47], which may result in defining an area with optimal characteristics for this species as an area of absence, could offer an incorrect outcome and decrease biological significance in the interpretation of the model [93,94].
In the MaxEnt configuration, 1000 was established as the maximum iteration parameter; the convergence limit was set at 0.00001 with a value of 0.001 for regularization [80]. The precision of the predictions was evaluated by using the area under the curve (AUC) [95]. AUC values below 0.7 were considered as poor, values between 0.7 and 0.9 were moderate and >0.9 were considered as high [96]. Hinge and threshold functions were disabled, leading to easily interpretable response curves and further adjustment to the ecological niche theory [97] (Supplementary Materials Figure S1). The program produced a set of maps in which each pixel represents a value between 0 and 1, the closer values to 1 indicate the greater suitability for the species, being thus interpreted as an index of habitat potentiality. It also allowed for the generation of response curves for the species in relation to the variables used, noting their suitability, the evaluation of optimal values, the tolerance intervals and the various variables' thresholds.

Presence Data
Presence data were collected through a huge bibliographical search combined with intensive fieldwork developed between 2011 and 2018 in Almería, Granada, Malaga, Murcia and Alicante provinces. Field data were geo-referenced using a GPS device (Garmin GPSMAP 60CX, 2 m error). Furthermore, digital sources were checked to gather distribution information on the species, such as GBIF [13], ANTHOS [98], and FAME project (Database of Threatened Flora of Andalusia) [99], and adding details about herbaria records (HUAL, GDA, GDAC, JAEN, MGC, and MUB). The bibliographical data, particularly those taken from the Internet, were carefully checked by QGIS [100] and using the latest aerial orthophoto, to only process reliable information. In order to avoid errors and remove duplicated occurrence records all the information was carefully filtered. A total of 1819 geo-referenced presence records of M. europaea were selected ( Figure 1). The dataset of presence records was analyzed by using a spatial autocorrelation analysis. the ecological niche theory [97] (Supplementary Materials Figure S1). The program produced a set of maps in which each pixel represents a value between 0 and 1, the closer values to 1 indicate the greater suitability for the species, being thus interpreted as an index of habitat potentiality. It also allowed for the generation of response curves for the species in relation to the variables used, noting their suitability, the evaluation of optimal values, the tolerance intervals and the various variables' thresholds.

Presence Data
Presence data were collected through a huge bibliographical search combined with intensive fieldwork developed between 2011 and 2018 in Almería, Granada, Malaga, Murcia and Alicante provinces. Field data were geo-referenced using a GPS device (Garmin GPSMAP 60CX, 2 m error). Furthermore, digital sources were checked to gather distribution information on the species, such as GBIF [13], ANTHOS [98], and FAME project (Database of Threatened Flora of Andalusia) [99], and adding details about herbaria records (HUAL, GDA, GDAC, JAEN, MGC, and MUB). The bibliographical data, particularly those taken from the Internet, were carefully checked by QGIS [100] and using the latest aerial orthophoto, to only process reliable information. In order to avoid errors and remove duplicated occurrence records all the information was carefully filtered. A total of 1819 geo-referenced presence records of M. europaea were selected (Figure 1). The dataset of presence records was analyzed by using a spatial autocorrelation analysis.  Table S1.

Environmental Variables
WorldClim [101] provides data on 19 bioclimatic variables considered as broadly significant from the biological point of view when modelling distribution areas [102][103][104]. Past bioclimatic data for the Mid-Holocene (ca. 6 ka AP), Last Glacial Maximum (LGM) (ca. 22 ka AP), and Interglacial Maximum (ca. 120 ka AP) were used in the present study. For the current distribution modelling, bioclimatic variables were generated from the  Table S1.

Environmental Variables
WorldClim [101] provides data on 19 bioclimatic variables considered as broadly significant from the biological point of view when modelling distribution areas [102][103][104]. Past bioclimatic data for the Mid-Holocene (ca. 6 ka AP), Last Glacial Maximum (LGM) (ca. 22 ka AP), and Interglacial Maximum (ca. 120 ka AP) were used in the present study. For the current distribution modelling, bioclimatic variables were generated from the climatic data period . The future estimation (2070) assumed the four Representative Concentration Pathways for emissions (RCPs), according to the fifth report from the IPCC [105]; RCP2.6 (a stringent pathway that requires that carbon dioxide (CO 2 ) emissions start declining by 2020 and reach zero by 2100, RCP4.5 (an intermediate scenario that requires emissions to start declining by approximately 2045 to reach roughly half of the levels of 2050 by 2100), RCP6.0 (a scenario where emissions peak around 2080 and then decline), RCP8.5 (it represents the basis for worst-case climate change scenarios, where emissions continue to rise throughout the 21st century). Bioclimatic data were downloaded from the Community Climate System Model (CCSM4) as the main reference for the distribution tests. This model simulates the global climate using four separate sub-models for the atmosphere, earth, oceans and sea ice [106]. These bioclimatic data had a 30-s resolution, except for the LGM, for which only data at 2.5 min were available. Each environmental variable map was adjusted to the Iberian Peninsula mask without Portugal.
To rule out the multicollinearity on the bioclimatic variables, a variance inflation factors analysis (VIF) was used. VIF analysis calculates variance-inflation and generalized variance-inflation factors for linear, generalized linear, and other statistical models able to discriminate and select the variables [107]. VIF analysis was performed with the R software [108][109][110]. According to this methodology VIF value must be under 5; thus, six bioclimatic variables were selected (BIO2, BIO3, BIO8, BIO9, BIO15, BIO19) to execute the modelling process (Table 1).

Potential Distribution of M. europaea
MaxEnt models' results for M. europaea distribution were adjusted by comparing its limits with the presence-absence of species that characterize different habitats together with [6,25,48]. Such species were Periploca angustifolia Labil. that coexist in the association Mayteno europaei-Periplocetum angustifoliae Rivas Goday and Esteve; Ziziphus lotus (L.) Lam, both species subsist in the association Gymnosporio europaei-Ziziphetum loti F. Casas, and mesophilic species, such as Buxus balearica Lam. and Cneorum tricoccon L., these three species coexist in the association Cneoro tricocci-Buxetum balearicae Rivas Goday and Rivas Mart. Chorological data for these plant species were obtained from GBIF and ANTHOS databases [13,98] (Figure S2). In addition, the current final map that describes the potential habitat of M. europaea and its decrease in south and southeast of Spain were designed as contiguous areas by combining presence records, species modelling, and photointerpretation of polygons performed in the national SIOSE project [111] as was done by Mendoza-Fernández et al. [47].

M. europaea SDMs in the Present
In the context of the present situation, the results for the M. europaea potential extend of occurrence (EOO) are expressed in two ways; values format by territory (see Table 2), and detailed distribution map (Figure 3).

M. europaea SDMs Towards the Future (2070)
Plots of the potential habitat projected into the future according to the representative concentration pathways are illustrated in Figure 4.

Distribution Dynamics
The model obtained for the LGM showed a significant dominance of the M. europaea habitat towards the southeast and the east of the Iberian Peninsula (Figure 2), exhibiting high suitability in areas nowadays submerged, and along the coastal plains of the Almería province. Probably, the populations were more stable in said area, acting as a genetic diversity reservoir in the glacial period. Médail and Diadema [112] identified the existence of such shelters as a crucial event in a context of global change. Furthermore, since the LGM possibly represented the most favorable scenario for this plant, the idea that M. europaea could not be as xerophyte an element as considered may be reinforced. During this period, weather conditions became more humid than in the Mid-Holocene and the Interglacial Maximum. In addition, in the case of the southern Iberian Peninsula, the temperature decrease was not as extreme as in other parts of Europe, favoring a climate characterized by relatively mild temperatures and intense rainfall [113]. On the contrary, the resulting model for both the Mid-Holocene and the Last Interglacial periods did not show such habitat suitability for the plant. While the results for both periods were quite similar, it should be noted that the Mid-Holocene scenario presented a medium suitability for this plant since the modelling process was very strict with the bioclimatic variables developed for this time interval. The outcome of this period subsequent to the LGM, where the niche model showed the greatest habitat amplitude for M. europaea in the southeast of Spain, may be understood as the consequence of a delay in the dynamics of the M. europaea populations from the LGM through the Mid-Holocene period, as previous studies in Sierra de Gádor suggest [23,114], which shows the gradual decline of M. europaea and other deciduous Quercus genus species that might have been more widely distributed during the LGM.
Otherwise, the results achieved by combining presence records, species modelling, and photo-interpretation of polygons demonstrated that M. europaea's current extent of occurrence (EOO) in Spain is approximately 166,721.2 ha, less than 47% of the suitable potential area (Table 2 and Figure 3). As shown in Figure 1, there is no presence of M. europaea in the east coast of the province of Almería or in a significant part of the Murcian coast. This is an important point to take into consideration when interpreting the modelling results, since this area is the driest in southern Spain, but where Z. lotus and P. angustifolia plant communities are able to develop and grow. Unifying the distribution of the aforementioned species in a single map helped to realize that M. europaea exhibits a less xerophytic behavior than Z. lotus and P. angustifolia, since the former does not coexist with the latter in the Murcian-Almeriensian bordering zone, thus proving M. europaea's more mesophilic character. The slight probability of presence recorded to the north and west of Spain may be caused by common temperature variables or the coincidence of the summer drought period, related to the savannoid character of this plant, since there is no evidence of any past or current presence in these areas. the aforementioned species in a single map helped to realize that M. europaea exhibits a less xerophytic behavior than Z. lotus and P. angustifolia, since the former does not coexist with the latter in the Murcian-Almeriensian bordering zone, thus proving M. europaea's more mesophilic character. The slight probability of presence recorded to the north and west of Spain may be caused by common temperature variables or the coincidence of the summer drought period, related to the savannoid character of this plant, since there is no evidence of any past or current presence in these areas.     In order to achieve higher accuracy in establishing the influence of bioclimatic scenarios, models of singular thermophilic and xerophilous species that coexist with M. europaea could be compared in equal space-time conditions [113]. Hence, several models from different communities may be obtained, and thus their spatial-temporal dynamics  In order to achieve higher accuracy in establishing the influence of bioclimatic scenarios, models of singular thermophilic and xerophilous species that coexist with M. europaea could be compared in equal space-time conditions [113]. Hence, several models from different communities may be obtained, and thus their spatial-temporal dynamics analyzed [86,87]. In addition, indirect gradient variables corresponding to the physical characteristics of the territory (elevation, orientation, slope, geology, etc.) could be included in the model (at least at present and future predictions), since they might show a good correlation with the patterns of species distribution by combining resource gradients and direct gradients [115].
Finally, the results of the future models with a projection for the year 2070 showed a slight variation in suitability among them, in which neither considerable habitats increase nor a significant decrease of it was observed. In the different scenarios, the zones exhibited small geographical variations, which slightly expand or reduce M. europaea's ecological niche (Figure 4). However, habitat suitability was not altered, remaining permanent, with intermediate values in the four models. The results reinforce the thesis of M. europaea's mesophilic behavior, and contrary to expectations, the global temperature increases projected in the different RCPs scenarios could predict the probable deterioration of the habitat. Nevertheless, a specific area was highlighted, the coastal plain in the south of the province of Almería (from sea level to 350 masl. approx.), where the four models corresponding to each RCPs scenario remained constant in terms of suitability, and offered the maximum values for all models. This fact makes this area a very probable one where potential M. europaea habitat is predicted in any of the 2070 scenarios. In addition, it supports the idea that this would be the most suitable area for the plant, in light of the four possible futures with an irradiative effect increase. Despite the fact that this area might be treated in the short term as one of the few habitats in the European continent available for this plant, it is currently a very altered territory, especially due to land-use changes resulting from the proliferation of tropical crops and intensive agriculture even in protected natural areas (an example of the deficient state of conservation in the Site of Community Importance (SCI) Artos de El Ejido (in Campo de Dalías, Almería, Spain) is shown in Figure 5), which has reached up to 90% of the potential area of this species [25,32,42,47].

Evolution of Distribution Patterns over Time
From a paleoecological point of view, the glaciation conditions could favor the extension of the potential niche for M. europaea in the Iberian southeast, corroborating that this plant does not present as xerophytic a behavior as Z. lotus or P. angustifolia, since the LGM model is the one that presents the highest suitability. This fact may be related to an increase in humidity in that period [116]. This supposition is consistent with two other types of evidence: a.
The presence of M. europaea growing and showing its maximum abundance dovetailing with that of the deciduous forests located in Sierra de Gádor during the Mid-Holocene according to palynological studies, b.
The link between M. europaea with C. tricoccon and B. balearica, which have mesothermophilic characteristics, capable of generating debate on the Maytenus-Periploca-Ziziphus trilogy, which was considered until now as a set of Ibero-African species with a hyperxerophytic conduct. This could explain the absence of M. europaea from Cabo de Gata to the Murcia region, being this area the most arid and with the smallest rainfall in the entire Spanish southeast, where conversely both Z. lotus and P. angustifolia are frequent.
In the near future, still affected by global change, the southern part of the province of Almería may become the most suitable area for M. europaea in the Iberian Peninsula. However, the natural area available for this plant s development might be limited by the proliferation of intensive crops, which already occupy more than 90% of the habitat in this zone. This might cause a serious conservation problem for this species and increase its extinction risk due to massive occupation, fragmentation and land-use changes.

Considerations for Conservation of This Plant Species in an Agricultural Matrix
Preservation of some semi-natural strongholds for the Spanish M. europaeea communities in the area known as Campo de Dalías (southern portion of Almería province) is fundamental for the conservation of the species. In addition, there are two additional strategies that may favor its preservation. The first one is related to the use of this species with the aim to establish hedges, vegetal fences [19] or small clusters that serve as a refuge for pollinators and predatory species of insects or other harmful pests for crops. In general, insects are considered one of the most effective animal groups in the pollination processes, thus their presence in natural systems is essential. Nevertheless, intensive agriculture and continuous landscape modification through land-use changes are associated with the decrease in the population of pollinators [117], and/or severe alterations in their population concentrations [118]. Moreover, the agricultural use of pesticides and broad-spectrum chemicals produces a decrease not only in pest insects but also in all the entomofauna associated with these systems, thus affecting the insects of less altered neighboring habitats [119]. The combination of Maytenus and Ziziphus bunches could considerably improve this strategy, and hence favor a more sustainable development model in an area where highly intensive agriculture production can be understood as an industrial system if inputs and residues are taken account of [25]. The second strategy of great interest would be the creation of green spaces, by means of peri-urban parks, where some of the remaining fragments of this community that cannot be conserved within a legally protected area, may be integrated.

Conclusions
We can conclude that although the Campo de Dalías area has suffered the greatest habitat loss in terms of the extent of occurrence and occupancy area, this is a key territory for the future of the species, as demonstrated by the distribution models generated throughout this research. Therefore, safeguarding the last remaining redoubts of the Spanish M. europaea communities in this zone is considered absolutely essential.
Moreover, it is revealed that the already deficient conservation status of M. europaea, a unique plant species in Europe, could worsen if the current global change drivers were to intensify, although the greatest threat factor for this species continues to be the loss of available habitat as a result of changes in land use. The limits of this research may be due to the enormous asymmetry in terms of spatial information of the taxon presence. This fact is most likely due to the deforestation of the coastal areas and the intense land changes that the south and southeast of Spain have suffered. Future lines of research could be aimed at observing how changes in climate are affecting the phenology of this species, designing ex situ conservation techniques and translocation of individuals, and at improving knowledge to determine if M. europea is in fact a relict "savannoid" paleotropical element rather than a "xerothermic scrub", and further detailing its distribution in Spain in order to guarantee the preservation of the highest number of populations and the maximum genetic diversity.