Recovery of Neglected Species with Cloud Water Micro Condense Capacity as a Response to Climate Change: The Case of Sclerophyllous Boxwoods of Buxus balearica Lam. in the Southern Spanish Mediterranean

Abstract

Understanding the survival needs of neglected or underutilized species (NUS) in agroforestry can offer valuable insights to mitigate climate change and biodiversity loss. This paper analyses the habitat heterogeneity of Buxus balearica populations (palaeorelict of the Mediterranean Basin) using a methodology that integrates four approaches: ecological profiling, multivariate analysis, and the clustering of populations according to environmental variables and suitability. The habitat analysis was conducted on 240 occurrence points, integrating open-source databases and germplasm collection field work. Results reveal that the distribution of B. balearica in the south of the Iberian Peninsula is mainly conditioned by thermal variations, the humidity regime, and microtopography, which makes it very vulnerable to current and future climate change scenarios. The analysis of habitat heterogeneity demonstrates its capacity to estimate genetic diversity, which provides a solid basis for future research and conservation actions. The rapid extinction process of these boxwoods is evident by comparing the suitability model under different climate scenarios. There is an urgent need to improve the current status of the species and associated landscapes, considering B. balearica as an underutilized species, and establishing a sound management plan to reinforce its populations and protect its natural habitat.

1. Introduction

The NUS (neglected and underutilized species) paradigm has become one of the most promising tools in the last two decades for policies to increase agrodiversity, sustainability, food autonomy, agricultural innovation, and adaptation to climate change [1]. The species and varieties that have the character of NUS, both neglected and forgotten, salvageable from ancient food and agricultural cultures, as well as underutilized, reduced in their use to small geographical or cultural areas, below their potential capacities, make up a block of forgotten genetic resources that are, in many cases, at a risk of genetic erosion or disappearance [2]. Their conceptual framework goes beyond the agri-food horizon to include hundreds of ornamentals, industrial, medicinal, and forestry species. Among the latter, we include those that, due to their ecosystemic role, can reduce water pollution, prevent soil erosion, reduce the CO₂ footprint, or improve aquifer recharge. For the last of these functions, important plant communities such as Laurel forests, Yungas, Guarani Atlantic forests, and Afromontane forests contribute to precipitating atmospheric water vapor [3,4,5,6,7]. Also, under Mediterranean climates and especially in their geographical environment, i.e., relict formations, species from the Late Tertiary vegetation (Messinian to Pliocene) survive, which retain a cryptoprecipitating character thanks to their evergreen, glabrous, microlauroid leaves, with oleoid, buxoid or myrtoid forms. They belong to genera such as Pistacia, Buxus, Rhamnus, Phyllirea, Coriaria, Maytenus, Osyris, Arbutus, Olea, or Quercus. Buxus balearica Lam. is one of these species [8].

It is, therefore, a residual species that continues to suffer a sharp reduction in its distribution area, frequency and abundance, a species that is disappearing along with other constituent elements of its thermophilic thickets (Maytenus senegalensis (Lam.) Exell., Osyris quadripartita Salzm. ex A.DC., Cneorum tricoccon L.), related to the lost laurel forests of the Mediterranean coastal mountains of the late Cenozoic. The profound impact of human presence on the environment of the riparian regions of the Mediterranean basin, forest fires and changes in the use and management of the forest environment, the disappearance of vectors related to their dispersal [9,10], and the accelerated climate change, generating the increase in summer water deficit, make the current existing protection measures insufficient. The areas declared as protected do not encompass all the small populations of boxwood in the Balearic Islands. While ex situ conservation measures such as seed banks aim to supplement these efforts by preserving boxwood accessions for the long term, the absence of reforestation or reintroduction plans for the species within its native habitats renders significant progress unattainable. B. balearica is not the only species of the genus threatened with extinction. This is also the case, for example, with B. vahlii Baill., a species endemic to Puerto Rico and Santa Cruz [11], or B. sinica (Rehder & E.H. Wilson) M. Cheng in China [12].

The leaf morphology of species such as B. balearica suggests the ability to condense moisture from fog and dew. Recent studies have demonstrated the influence of the Atlantic fog regime on the distribution of populations of Argania spinosa (L.) Skeels, with similar leaf morphology, in their natural range (mainly in Morocco), shows remarkable habitat heterogeneity [13,14]. A similar phenomenon is observed in B. balearica, according to the data and experience we have acquired in the last 30 years after monitoring its populations and studying the habitat and ex situ conservation of this species in the south of the Iberian territory [15,16,17]. The relict character and its presence in two different climatic floors, the thermo- and lower meso-Mediterranean, have been already confirmed by its phytosociology [18,19]. There is also evidence of genetically different populations at the molecular level associated with this phytogeographical heterogeneity [15,20].

B. balearica still occurs in the Western Mediterranean, in Sardinia, Balearic Islands (in Cabrera and Mallorca and in the Sierra de Tramontana), and very locally in N. Morocco (Rif, Middle Atlas, Saharan Atlas, and High Atlas). We also find possible historical references in other Tyrrhenian territories, as the commentary of Theophrastus (3rd century BC) refers to various localities and uses of boxwood, speaks of the very large and splendid boxwoods of Cirno (Corsica), and indicates how those of this island are more robust than those of any other place, which leads us to think that it may be B. balearica.

In the Iberian Peninsula, B. balearica is found exclusively in Andalusia (the provinces of Malaga, Granada, and Almeria). It occupies some coastal populations and other more inland ones until 1.900 m.a.s.l., between the mountains of Cázulas, Los Guajares, Tejeda, Almijara, and Alhama. It is also found in the Sierra de Gádor [17,21,22,23,24]. The area is subject to a climate that varies from the semi-arid Mediterranean (below 500 m.a.s.l.) to the sub-humid Mediterranean (above 600 and 700 m). Average monthly temperatures range from 11 °C to 26 °C and can reach up to 40 °C during the summer. Rainfall ranges from 450 to 500 mm in coastal areas and over 1000 mm in some mountain areas.

The Mediterranean influence is characterized by the interaction of warm and humid sea air with the land surface, leading to the condensation of water vapor and the formation of coastal fog. The region’s abrupt topography further shapes this meteorological phenomenon, acting as a natural barrier that compels the moist air from the Strait of Gibraltar and the Alboran Sea to ascend, cool, and condense, giving rise to orographic fogs [25,26]. Also, the cloudy fronts that reach these coastal mountains from the peninsular interior, pushed by winds of the west, create when descending after overcoming summits and water dividers, another type of orographic fog that generates humidity at the highest altitude levels.

The present study aims to identify and corroborate the presence of B. balearica boxwoods in the two types of habitats mentioned in the Iberian area studied, one sheltered at high altitudes and the other one that is more thermophilic and coastal. Furthermore, we aim to assess the significance of advection fogs in facilitating the survival of this species, particularly in coastal environments.

Advection fogs, a meteorological phenomenon of interest, form when the prevailing warm and dry westerly winds from the interior, locally referred to as “terral”, cease and push warm surface waters towards the Alboran Sea. This, coupled with the subsequent thermal inversion, condenses substantial humidity masses when the moist air comes into contact with a colder surface. These condensed fogs are transported inland, affecting coastal regions [25,26]. Remarkably, these fogs have been documented for millennia, with Phoenician colonies inhabiting this territory from the 8th to 3rd centuries BC. The Phoenicians referred to these fogs as “taró”, a term that has persisted up to the present day. These ancient mariners made offerings to the goddess Malak (Malaga, formerly Phoenician Malaka, lies within the sphere of influence of taró) in hopes of securing success in coastal fishing—a tradition that continues to hold economic and cultural significance in the region.

This article undertakes a comprehensive statistical and probabilistic analysis of the environmental factors that define the Andalusian region. Our goal is to present a thorough examination of the distribution patterns of B. balearica. Building on the suitability study previously conducted [17], we have integrated new climate change scenarios to account for periods of temperature fluctuations that would result in a greater water stress for the species. This information will inform the development of an ex situ conservation program that takes into account the assessed genetic diversity.

2. Materials and Methods

2.1. Case Study Species and Its Ex Situ Conservation Program

The fieldwork consisted of collecting thirty-seven new locations of B. balearica that were included in a germplasm collection programme for ex situ conservation and analyzing the limits of the areas currently protected by regional and national legislation in Andalusia and Spain. The thirty-seven accessions are preserved in a seed bank using cold storage techniques. These accessions were collected over the last ten years and based on a variability criterion by selecting individuals from as many types of habitats as possible [27,28].

Since the presence of fog significantly increases atmospheric saturation, relative humidity data have been obtained from the ASOS (Automated Surface Observing System) network, using the station closest to the study area that is located in Málaga airport at latitude 36.67 N and longitude 4.49 W (LEMG) (https://mesonet.agron.iastate.edu (accessed on 17 October 2021)). Although other stations within the study area are part of the national observation networks operated by the State Meteorological Agency (AEMET), we exclusively relied on the ASOS station with 30 min data availability. This specific time resolution is essential to discern anomalies in records potentially linked to abrupt shifts in humidity that could result from fog events.

Source Data

Habitat analysis of B. balearica has been carried out using 240 occurrence points in the south of the Iberian Peninsula (provinces of Malaga, Granada, and Almeria) obtained from the GBIF and REDIAM databases, and mainly the accessions collected for their subsequent long-term ex situ conservation in the Andalusian Plant Germplasm Bank.

The information layers have been downloaded from the environmental information portal in Andalusia (https://portalrediam.cica.es/descargas, accessed on 30 July 2021), from which 20 environmental variables of temperature, precipitation, topography, and insolation have been obtained (

Table 1. Climatic and environmental variables collected for this ecogeographic study.

Code	Variable	Description
TMC	Mean temperature of the warmest month (°C)	Mean value in August * Maximum value in August *
TMMC	Maximum temperature of the warmest month (°C)
TMF	Mean temperature of the coldest month (°C)	Mean value in January * Minimum value in January *
TMMF	Minimum temperature of the coldest month (°C)
PT01	Mean precipitation of the first quarter (mm)	January–February–March *
PT02	Mean precipitation of the second quarter (mm)	April–May–June *
PT03	Mean precipitation of the third quarter (mm)	July–August–September *
PT04	Mean precipitation of the fourth quarter (mm)	October–November–December *
MDT	Elevation from the Digital Elevation Model (m)	DEM of Andalusia
ORI	Aspect (degree)	Calculated from DEM
SLP	Slope (%)
DIST	Euclidean distance to the drainage network (m)	Calculated from water bodies
MDI1	Intermediate insolation model	Between December–March
MDI2	Insolation model during spring equinox	On 21 March
MDI3	Intermediate insolation model	Between March–June
MDI4	Insolation model during summer solstice	On 21 June
MDI5	Intermediate insolation model	Between June–September
MDI6	Insolation model during autumn equinox	On 21 September
MDI7	Intermediate insolation model	Between September–December
MDI8	Insolation model during winter solstice	On 21 December

* The period from 1971 to 2000.

).

The variable values for each occurrence point were extracted using the QGIS point sampling tool plugin that combines the information layers and the points of presence in the Coordinate Reference System: ETRS89 UTM 30 N, with a pixel size of 75 m for each variable.

To perform the three quantitative analyses, we standardized the input data using the z-score formula, ((Xi − X)/s), where Xi represents a specific data point, X is the sample mean, and s represents the sample’s standard deviation.

2.2. Numerical Ecology Approach

2.2.1. Ecological Profiles

The ecological profiles allow us to check the ecological requirements of the study species by analyzing the frequency distribution of the environmental data layers. The comparison of the maximum and minimum values for the entire range allows for the identification of ecological behaviors and preferences. Each variable acts as a descriptor, and for its representation, the generic function of R has been used, hist() [29].

2.2.2. Multivariate Analysis

The correlation analysis between the environmental variables has been achieved by applying the principal component analysis (PCA), which determines the variance explained by these few independent main axes and the contribution of the variables in defining these axes. The variables that provide a more significant variance in the construction of the axes have been used in the clustering analysis.

A hierarchical fusion method has been selected to find groups of individuals with similar demands and to differentiate populations. Clusters were built using the Ward.D2 method that was executed with the hclust() function of the STATS package of R version 1.3.1093.

2.2.3. Suitability Models

Climate change models used to predict the potential species distribution were the Local Climate Change Scenarios of Andalusia (ELCCA) updated to the 4th Report of the Intergovernmental Panel on Climate Change (IPCC). These models are available in the catalogue of data for Andalusia in three climatic periods, i.e., 2011–2040, 2041–2070, and 2071–2099, calculated from the combination of the general circulation models of the atmosphere (CGCM2 and ECHAM4) and different emission scenarios (A2 and B2). These data were downloaded from the environmental information portal in Andalusia. The climate model used was the BCM2 in the sra2 scenario corresponding to average temperature values of 16.7 °C, 17.3 °C and 18.9 °C in the periods 2011–2040, 2041–2070, and 2071–2099, respectively [30].

The three potential scenarios have been modelled with the maximum entropy algorithm of MaxEnt software version 3.4.3, November 2020 [31,32]. The occurrence points provide the values of each variable in both the current and future scenarios described above. In this way, it is possible to predict the regions with the highest probability of suitability that meet the maximum entropy values. The suitability scale is represented between 0 (lowest suitability) and 1 (highest suitability), and the area under the ROC curve (receiver of characteristics) shows the discriminative capacity of the calculated model and presents a tighter distribution if greater than 0.75 (poor 0.75 < AUC > 0.95 excellent). The model configuration parameters have been calculated using the ENMevaluate function of the ENMeval package of R (v. 2.0.3) [33].

3. Results

The process of collecting accessions for ex situ conservation has enabled the incorporation of the accessions listed in

Table 2. Accessions of Buxus balearica (Andalusian Seed Bank (BG), Spain).

BG Accession Number	Collection Locality	Province	Type of Collection
14779	Nerja	Málaga	Base
15596	Frigiliana	Málaga	Base
15636 *	Nerja	Málaga	Base
16208	Nerja	Málaga	Base
16662	Vélez de Benaudalla	Granada	Base
17667	Vélez de Benaudalla	Granada	Base
17779	Almijara	Málaga	Base
18579	Almuñécar	Granada	Base
19213	Motril	Granada	Base
19214	Lobres	Granada	Base
19215	Cómpeta	Málaga	Base
19216	Almijara	Málaga	Base
19217	Cómpeta	Málaga	Base
19218	Cómpeta	Málaga	Base
19219	Nerja	Málaga	Base
19220	Cómpeta	Málaga	Base
19221	Cómpeta	Málaga	Base
19222	Cómpeta	Málaga	Base
19223	Cómpeta	Málaga	Base
19224	Nerja	Málaga	Base
19225	Ragor	Almería	Base
19226	Cázulas	Granada	Base
19227	Otívar	Granada	Base
19238	Arroyo de la Miel	Málaga	Base
19520	Güájares (Los)	Granada	Base
20593	Vélez de Benaudalla	Granada	Base
22061	Vélez de Benaudalla	Granada	Base
50949	Guadalfeo (ex situ)	Granada	Active
51043	Guadalfeo (ex situ)	Granada	Active
51439	Guadalfeo (ex situ)	Granada	Active

Base = −20 °C storage condition. Active = −5 °C storage condition. * Accession 15636 is stored in Active and Base collections.

into the collections of the Bank of Andalusian Plant Germplasm (Junta de Andalucía) collections.

From 2004 to 2021, during the warm months of June, July, and August (refer to

Table 3. Historical (2004–2021) summary of records of total hours in June, July, and August with relative humidity obtained from the LEMG station (ASOS network 36.67 N, 4.49 W) above 70% (comfort point) and 90% (condensation point).

	RelH > 70%						RelH > 90%
Year	Records	June	July	August	Total Hours	%	Records	June	July	August	Total Hours	%
2004	1586	434	676	476	793.0	37	65	14	9	42	32.5	2
2005	1025	279	285	461	512.5	24	4	3	0	1	2.0	0
2006	1175	350	317	508	587.5	27	3	0	0	3	1.5	0
2007	1164	276	420	468	582.0	27	14	1	4	9	7.0	0
2008	1316	361	430	525	658.0	30	31	1	24	6	15.5	1
2009	1141	352	421	368	570.5	26	21	19	0	2	10.5	0
2010	1615	369	545	701	807.5	37	62	38	18	6	31.0	1
2011	1786	656	480	650	893.0	41	70	30	5	35	35.0	2
2012	1561	433	489	639	780.5	36	233	91	38	104	116.5	5
2013	1481	372	431	678	740.5	34	40	17	14	9	20.0	1
2014	1838	460	620	758	919.0	43	443	74	177	192	221.5	10
2015	1173	215	580	378	586.5	27	56	16	40	0	28.0	1
2016	1763	506	521	736	881.5	41	37	14	16	7	18.5	1
2017	1779	471	574	734	889.5	41	58	9	32	17	29.0	1
2018	1096	271	312	513	548.0	25	59	39	8	12	29.5	1
2019	1449	404	426	619	724.5	34	69	37	3	29	34.5	2
2020	1583	305	744	534	791.5	37	89	19	29	41	44.5	2
2021	1824	483	581	760	912.0	42	126	36	18	72	63.0	3
Total	26,355	6997	8852	10,506	13,177.5		1480	458	435	587	740.0

and Figure 1), a total of 13,177 accumulated hours of relative humidity above 70% were recorded. Additionally, there were only 740 days with humidity exceeding 90%, with August showing the highest humidity in both cases. This observation may align with the increased fog frequency in the Alboran Sea during the warm period (June–August). It is asserted that during periods of heightened water stress, there is a noteworthy contribution of humidity in the area.

Observations reveal cyclical fluctuations occurring approximately every 7 years. Margalef [34,35] also references these climatic cycles in connection to the occurrence of colder winters and their impact on fishing yields along the Mediterranean coasts. This pattern could potentially be associated with advection fog events.

As can be seen in the temperature profiles (Figure 2), in the areas with a range of 4–16 °C, the species prefers lower temperatures, between 8–10 °C, and dropping to 4 °C during the coldest month. Regarding higher temperatures, there are no significant differences when compared to the rest of the territory.

The bimodal distributions of the frequencies related to temperatures (Figure 2), precipitation (Figure 3), altitude (Figure 4), and insolation models (Figure 5) show the existence of two populations.

The precipitation graphs by quarters (Figure 3) reveal a distinction between winter and autumn, suggesting that the existence of drier localities is possibly related to Almería. Notably, there is a preference for areas with higher precipitation. During winter, the species aims to accumulate the greatest amount of water (200–300 mm). Remarkably, in the summer, the species is found in areas with lower rainfall, indicating a resistance to drought or reliance on alternative sources of humidity beyond precipitation.

As shown in Figure 4, altitude presents a Gaussian distribution, tending to be bi-nomial, with a decrease in frequency between 300 and 500 m. This suggests the existence of populations favoring coastal areas and others with an optimum of around 800 m, disappearing above 1200 m. Orientation and slope show a slight preference, or a better preservation, with southern orientations.

The distance from the axis of the drainage network indicates that the plant tends to thrive in areas near ravines, streams, and rivers. This preference is likely due to lower incident radiation, reduced evapotranspiration, and enhanced defense against forest fires.

During equinoxes, the distribution exhibits a greater variance (Figure 5), allowing for a better diagnosis of heterogeneity. However, a stronger correlation with the distribution of localities is observed during solstices. Specifically, boxwood localities prefer areas with fewer hours of sun in winter (minimum sun elevation) and more sun in summer (maximum sun elevation). The MDI4 and MDI8 models best explain the distribution of populations.

There is a distinct behavior in the variables MDI12 and MDI22 compared to MDI0, evident in the factorial study. The ecological profile reveals that during days close to the winter solstice, numerous localities do not receive direct sunlight. Conversely, the profile during the summer solstice should feature many localities with a high number of sunshine hours, aligning with the equinox profiles. The October–November profile is intermediate.

These results correspond to those observed in the factorial study and the cluster of the species’ locality distribution. Specifically, the temperature frequency distributions show a clustering of localities in a thermophilic coastal sub-population and in another of a montane character, occupying the meso-Mediterranean level. The species does not seem to tolerate frost (or at least does not respond well to it due to competition with other species), disappearing below the mean temperatures of the coldest month below 7 °C and mean minimum temperatures of the coldest month below 3 °C.

According to PCA analysis, component 1 explains 38% and comprises the variables of the insolation models, with the spring and autumn equinoxes contributing the greatest weight (Figure 6). In contrast, the variables MDI4 and MDI8 contribute 0.26 and 0.32, respectively, to the same component (

Table 4. Loadings and the accumulated variance explained with principal components.

Code	Variable	Principal Component (PC)
		1	2	3
TMC	Mean temperature of the warmest month (°C)	0.12	−0.34	0.07
TMMC	Maximum temperature of the warmest month (°C)	0.05	−0.22	0.40
TMF	Mean temperature of the coldest month (°C)	0.13	−0.29	0.35
TMMF	Minimum temperature of the coldest month (°C)	0.13	−0.33	0.05
PT01	Mean precipitation of the first quarter (mm)	−0.09	0.34	0.26
PT02	Mean precipitation of the second quarter (mm)	−0.09	0.34	0.20
PT03	Mean precipitation of the third quarter (mm)	−0.12	0.28	0.36
PT04	Mean precipitation of the fourth quarter (mm)	−0.09	0.32	0.37
MDT	Elevation from the Digital Elevation Model (m)	−0.12	0.31	−0.37
ORI	Aspect (degree)	0.24	0.01	−0.04
SLP	Slope (%)	−0.18	−0.10	−0.21
DIST	Euclidean distance to the drainage network (m)	−0.01	0.16	−0.33
MDI1	Intermediate insolation model	0.32	0.08	−0.09
MDI2	Insolation model during spring equinox	0.34	0.11	−0.04
MDI3	Intermediate insolation model	0.32	0.14	0.06
MDI4	Insolation model during summer solstice	0.26	0.15	0.12
MDI5	Intermediate insolation model	0.32	0.14	0.06
MDI6	Insolation model during autumn equinox	0.34	0.11	−0.04
MDI7	Intermediate insolation model	0.32	0.08	−0.09
MDI8	Insolation model during winter solstice	0.30	0.05	−0.10
Standard deviation		2.76	2.32	1.57
Proportion of Variance		0.38	0.27	0.12
Cumulative Proportion		0.38	0.65	0.77

). Component 2 groups in one direction the temperatures of the warmest month and the minimum temperatures of the coldest month and in another the precipitation from January to June with an accumulated variance of 27%. Component 3 only contributes 12% but absorbs the topographic variables (MDT and Dist), the precipitation from July to December (PT03 and PT04), the maximum temperatures of the warmest month (TMMC), and the mean of the coldest month (TMF).

The first three axes together account for 77%. Axis 1 is highly correlated with the degree of insolation (and shielding in the opposite direction) that the vegetation receives, especially with respect to the number of hours of sunshine at the summer and winter equinoxes. Axis 2 is mainly configured according to the extreme temperatures recorded (minimum in the coldest months and maximum in the hottest months). The configuration of axis 3 is dominated by precipitation in the second half of the year (August onwards) together with altitude.

The 12 variables that contribute the most variance (>|0.33|) to the three principal components are TMC, TMMC, TMF, TMMF, PT01, PT02, PT03, PT04, MDT, DIST, MDI2, and MDI6 (Table 4). However, those with a correlation higher than |0.75| (TMF, TMMF, MDI2, and MDI6) have been removed. Due to the contribution of the insolation models during the two solstices (MDI4 and MDI8) and given that there is no correlation between these variables, the eight chosen variables are TMC, TMMC, PT01, PT03, MDT, DIST, MDI4, and MDI8 (Figure 6).

The variables included in the cluster analysis were PT01, PT03, TMC, TMMC, MDT, DIST, MDI4, and MDI8 (Figure 6). There are two sets of well-differentiated localities, those of the lower meso-Mediterranean and upper thermo-Mediterranean levels that are very thermal and free of frost with average temperatures of the coldest month above 11 °C and those of the upper thermo-Mediterranean meso-Mediterranean with average temperatures that drop to 7 °C.

According to Figure 7, within the coastal and thermal sub-population, there are five different nuclei: two more eastern ones (Rágol and Guadalfeo–Molvízar), another central one (Cantarriján), another western one (Nerja and lower Frigiliana), and finally the westernmost nucleus (Cómpeta). Within the mountain sub-population, we can distinguish the nucleus of Lentejí and Otívar, Competa and upper Nerja, and Los Guájares.

In these results (Figure 8), the Balearic populations were immediately differentiated from those of the Iberian Peninsula, and within these, the eastern population of the Sierra de Gádor (Almeria) was separated from the Malacitano-Granadina populations, where again the coastal populations were separated from the inland mountains, and in turn, these were separated in a western block (Sierra de Almijara: Nerja-Frigiliana) from the eastern populations of the Sierras de Cázulas and Guájares (Guadalfeo-Molvízar).

Different models of potential presence calculated for B. balearica using Maxent in the periods 2011–2040, 2041–2070, and 2071–2099 (Figure 9) were obtained with a reliability of 0.98, according to the analysis of the AUC values of each one. The algorithm effectively determines habitats, particularly under different climate change scenarios.

In the 1961–2000 scenario, the potential distribution extends west of the current central core through the Malacitan mountains. The 2011–2040 scenario shows a slight variation, despite demographic regression in current populations. The 2040–2070 scenario reduces the density of presence in its central core in Málaga–Granada, and its potentiality decreases more drastically in the western potential core. In the 2070–2100 scenario, B. balearica landscapes disappear completely from their western and eastern potential cores, while their presence in the current central core becomes anecdotal. The species and its landscapes disappeared almost entirely from the Iberian Peninsula.

4. Discussion

The study of a species’ habitat was, for a long time, the traditional method of assessing genetic diversity in forest sciences until the advent of molecular techniques. The difficulty of conventional genetic methods and the longer life cycle of most forest species justified this form of assessment. Recent studies [26,36] with Argania spinosa have already demonstrated the enormous usefulness of this approach, which is again evidenced here with B. balearica. Our strategy simultaneously and synergistically uses several numerical methods, from the classical study of ecological profiles (statistical analysis) and clustering and principal component methods (multivariate analysis), together with predictive models that include different climate change scenarios (Maxent), allowing us to compare and correlate results with those obtained from the study of genetic diversity at the molecular level [16,20,37].

The frequency distribution of the reticules in which the species appears does not correlate with the frequency of precipitation in all the reticules studied, except in the summer quarter, where its presence seems related to a low precipitation (between 20 and 40 mm). These localities receive an annual rainfall of >500 mm, with about 110 mm in the second quarter. Therefore, the species seems to require a certain reserve of rainwater in the soil, to which a quantity not yet estimated should be added, coming from leaf condensation. Given the climatic variation observed in the region over the last few decades, this water supply may be becoming increasingly scarce.

The habitat of B. balearica in the Southern Spanish Mediterranean suggests the existence of two very different sets of localities: those of the lower thermo-Mediterranean floor, very thermal, frost-free, with average temperatures in the coldest month above 11 °C, and those of the upper thermo- and meso-Mediterranean with temperatures that even drop to 7 °C. Indeed, within the coastal and thermal sub-population, an eastern nucleus (Guadalfeo-Molvízar) can be clearly distinguished from two others, a central one (Cantarriján and Arroyo de la Miel) and a more westerly one (lower Nerja). Within the mountain sub-population, several nuclei can be distinguished: (a) the nucleus of Otívar, probably characterized by the high slopes (the boxwoods are on vertical walls); (b) the upper eastern nucleus of Guájares-Cázulas; (c) the central nucleus of the upper valley of the river Chíllar; (d) the western nucleus of the valleys above Cómpeta, including that of Cázulas; and (e) the small nucleus of Canillas de Aceituno (Figure 8).

Furthermore, it is corroborated that the species only occurs on limestone soils with a pH > 7. In the provinces of Málaga and Granada, they are always found on kakiritized marmoreal limestones; and in the Almería locality of Rágol, they are also found on limestones, although not marmoreal.

Similarly, due to the low dependence on topographic variables, it can be deduced that the species tends or is tending to take refuge in those areas where it has a low-temperature oscillation and optimum humidity without following any strictly altitudinal or latitudinal pattern, as is the more general case of other species.

It remains to be demonstrated whether the habitat diversity is only due to the disappearance of boxwood in intermediate situations, between 300 and 500 m altitude, due to external aspects (greater occupation or human pressure), or whether these two phenotypically and genotypically differentiated populations exist.

The resulting clusters coincide with those obtained by other investigations at the molecular level using RAPDs [16,20], estimating practically the same grouping of populations. The easternmost population in Almería (Rágol) is differentiated from that of the Malacitano-Granadino group, and within these, a coastal nucleus was distinguished from a montane (meso-Mediterranean) nucleus. Within the latter, there were also differences between the eastern localities (Lentejí, Otívar, Guadalfeo, and Molvízar) and the western ones (Cómpeta). The genetic diversity observed by other authors [16,20] was greater at the intra-population level than at the inter-population level (because of the cross-population nature of the species). However, a significant genetic difference was observed between the populations of the meso-Mediterranean group.

Concerning the potential distribution of B. balearica in future climate change scenarios, the habitat of B. balearica would not only be reducing and fragmenting, but also it would be lost in the strict sense of the term. The climatic conditions that allowed the survival of a forest made up of myrtoid, ramnoid, or buxoid leaf species, which are coriaceous, glabrous, non-spiny, and capable of achieving cloud water interception by condensing the microclimatic humidity of the soft mists and sea breezes that occur even in the middle of summer, would disappear with them. This species was not and is not the only species with this function, as others within the same community, such as Rhamnus alaternus L., Pistacia lentiscus L., Coriaria myrtifolia L., or Quercus coccifera L., are species with leaves or leaflets like those of Buxus [8,18,36].

The application of the model in several consecutive future climate scenarios shows a process of progressive and rapid extinction of the boxwoods of B. balearica in the south of the Iberian Peninsula. The area of its potential presence under the climatic conditions of the first quarter of the 21st century is drastically reduced in only twenty more years, becoming almost completely extinct in the eastern half of its current area. By 2070, their extinction would be already profoundly affecting the western half as well, and by the end of the century, their extinction would be total. Unfortunately, climate change predictions are coming true faster than forecasts based on probability calculations [38]. As Schultes [39] said, in relation to the loss of traditional knowledge associated with biodiversity and the need for urgent action, also with regard to the in situ and ex situ conservation of B. balearica as with many other NUS species, it will soon be too late.

All these boxwoods constitute cultural landscapes formed by B. balearica, and many of the species mentioned above are phytosociologically associated with it. These are NUS (neglected or underutilized species) of a forestry nature that, in the past, were used not only for their firewood, wood, mastic, and other resins but also for a wide range of aromatic, medicinal, pickling, dyeing, melliferous, and condimentary uses. In addition to these uses, the important ecosystemic function is represented by its crypto-precipitating capacity, condensing fog (advection and orographic), which represents a capacity to recharge aquifers that has unfortunately largely disappeared today.

As a result of the habitat diversity study that was carried out, 65 doses of Balearic boxwood seeds from 20 different localities in the provinces of Granada, Malaga, and Ameria have been collected and conserved in the Andalusian Plant Germplasm Bank (Table 2). Ex situ cultivation has already generated new accessions in the Royal Botanical Garden of Cordoba. This set constitutes a genetic reserve that ensures its conservation and possible future in situ actions using reintroduction techniques. The study could also seriously correct the limits of the current area protected by the declared Natural Park “Sierras de Tejeda, Almijara, and Alhama” since many of the localities studied and collected for ex situ conservation are outside these limits. This makes ex situ and in situ conservation measures essential to ensure the future of these boxwoods and other rare and threatened species and plant communities in the Mediterranean region [40,41]. This involves implementing a comprehensive plan that includes expanding the boundaries of protected natural areas, choosing specific locations for population reinforcement or reintroduction, and, crucially, fostering a nuanced understanding of these ecosystems’ value. Not only do they serve as cultural landscapes, but also, they play a pivotal role in cloud water interception and water management within the Andalusian Mountain range.

5. Conclusions

This comprehensive study of Buxus balearica, a species integral to the Mediterranean’s highly fragmented cultural landscape, underscores its status as a neglected and underutilized species (NUS) with significant forest interest. Facing threats of extinction due to anthropogenic pressures, marginalization, and climate change, B. balearica’s conservation is crucial. The ex situ conservation efforts have been successful, with 31 new accessions from diverse localities like Nerja, Cómpeta, Cantarriján, and Guadalfeo, including a novel entry from Rágol, Almeria.

The study’s habitat heterogeneity analysis, building upon previous genetic diversity estimates, reveals critical insights into the environmental needs of B. balearica. It highlights the species’ distinct distribution in the Iberian Peninsula, separating into two sets of localities: those from the lower thermo-Mediterranean level influenced by summer water compensation through advection fogs and those from the high-Mediterranean and meso-Mediterranean levels, dependent on orographic humidity at higher altitudes. This distinction is vital in understanding the species’ adaptation and survival strategies.

Significantly, the study shows the species’ sensitivity to summer water stress and its reliance on microclimatic phenomena, highlighting its vulnerability to the ongoing climate change in the region. The detailed ecological profiles, including thermometric and pluviometric variables, provide a nuanced understanding of the species’ preferred environmental conditions, showing a preference for areas with higher precipitation and lower temperatures.

The findings also reveal a bimodal distribution in the species’ frequency related to temperature, precipitation, altitude, and insolation models, suggesting the existence of two distinct populations. This is further elucidated through the analysis of principal components, which sheds light on the interplay between various ecological factors like insolation, temperature extremes, and precipitation patterns in shaping the species’ distribution.

The potential distribution models for B. balearica, developed using MaxEnt ecological niche modeling (software version 3.4.3, November 2020), indicate a significant shift in the species’ habitat under future climate scenarios. The predictive models, with high reliability, show a worrying trend of habitat contraction and loss, particularly in the species’ western and eastern potential cores, with a marked reduction in the central Málaga–Granada area.

In summary, this study not only contributes to the understanding of B. balearica’s ecological requirements and its response to climatic changes but also highlights the urgency for targeted conservation strategies to preserve this valuable Mediterranean species.