Vanilla distribution modeling for conservation and sustainable cultivation in a joint land sparing/sharing concept

. Vanilla, an expensive but popular spice used in many industries, faces problems related to its supply. Some of these problems are due to the fact that vanilla cultivation is based on clonal material of a single species ( Vanilla planifolia ) and is dominated by just a few countries located outside the native growing areas of aromatic vanilla species, which is the neotropics. Despite the economic importance of this crop, relatively little attention has been paid to its wild relatives, in particular with respect to their biology, ecology, and potential use. We hypothesized that species distribution models (SDMs) can identify suitable areas for both the conservation and cultivation of vanilla crop wild relatives (CWRs), following a joint land sparing/land sharing (SPASHA) approach, thus offering alternative sourcing areas and production methods. This is the ﬁ rst study that explored the use of ensemble SDMs to provide applicable land use maps related to the conservation and sustainable cultivation of wild vanilla species in Costa Rica, contributing to a solution for the problems related to current vanilla production systems. We focused on four aromatic vanilla CWRs, native to Costa Rica, to make land use policy recommendations for this country, and more speci ﬁ cally for the biological corridor Osa and its surroundings within our study region (cid:1) Area de Conservaci (cid:1) on Osa (ACOSA). The resulting distribution maps, with a mean AUC of 0.89, re ﬂ ected their current potential distribution (ranging from unsuitable to suitable) in Costa Rica. Combining them with recent land use and conservation area maps of our study region, we de ﬁ ned (1) areas for vanilla conservation and (2) areas for sustainable vanilla cultivation within agroforestry systems. These land use recommendations can now be integrated within the National Bio-Corridor Program (PNCB) that aims at making biological corridors more productive by proposing alternative income generation for local communities living within these areas. Our approach can be applied to identify priority areas for implementing the SPASHA approach on other vanilla CWRs and in more regions across its native growing ranges, given the availability of land use maps and enough occurrence records to build accurate SDMs.


INTRODUCTION
Vanilla is an expensive but crucial ingredient for many industries worldwide (Ramachandra Rao and Ravishankar 2002, Havkin-Frenkel and Belanger 2010, Cameron 2011. It is extracted from the fruits, better known as beans or pods, of orchid vines belonging to the genus Vanilla Plumier ex. Mill (1754) Dressler 2010, Cameron 2011). At present, the main production countries are Madagascar (3227 tons/yr), Indonesia (2402 tons/yr), and China (662 tons/yr; Faostat 2019), all three located outside the native growing range of aromatic vanilla species, which is the neotropics (Purseglove 1973, Bouetard 2010, Soto Arenas and Dressler 2010. Vanilla cultivation in the aforementioned countries is therefore limited to the introduced species Vanilla planifolia (Fouch e andJouve 1999, Lubinsky 2006). Due to the fact that vanilla orchids are propagated in a vegetative way, the genetic diversity of V. planifolia crop populations is relatively low, negatively affecting production resilience (Besse 2004, Bory 2010, Borbolla P erez 2016. Furthermore, the lack of natural pollinators in the above-mentioned countries forces vanilla producers to rely on hand pollination, a very labor-intensive and delicate activity (Villanueva-Viramontes 2017). As hand pollination results in higher pollination rates and thus yields compared to natural pollination, this technique has been applied in all commercial vanilla plantations worldwide, including those in the neotropics (Cameron 2011). Despite its high efficiency, it potentially leads to over-pollination, which stresses the plants and makes them less tolerant to abiotic and biotic stresses, such as droughts, temperature changes, and diseases caused by Fusarium, Phytophthora, and Glomerella and by the cymbidium mosaic virus (Havkin-Frenkel and Belanger 2010, Cameron 2011, Ramos-Castell a et al. 2016, Villanueva-Viramontes 2017. This downward spiral is exacerbated by the lack of genetic diversity of the vanilla orchids in plantations. The concentration of vanilla production in just a few countries also leads to a lack of production system's resilience against extreme weather events, pests and diseases, political instability, theft, monopoly, and processing errors (Havkin-Frenkel and Belanger 2010). Notwithstanding the economic importance of this crop and the problems related to current vanilla production systems, other aromatic vanilla species besides V. planifolia have received relatively little attention from the scientific community, in particular with respect to their biology, ecology, and potential use. More research is urgently needed to determine how we can preserve these wild species while also exploring the possibilities to cultivate them in sustainable production systems and improve current systems based on the single species V. planifolia (Ramos-Castell a et al. 2016).
The lowland rainforests of the neotropics are home to several aromatic vanilla crop wild relatives (CWRs), which are wild species related to the vanilla crop (Bory 2010, Flanagan andMosquera-Espinosa 2016). These CWRs, which are naturally pollinated (Lubinsky 2006, Pansarin and Pansarin 2014, Anjos et al. 2016, Pansarin and Miranda 2016, represent valuable sources of useful traits to improve current vanilla cultivation systems in terms of climate change adaptation, disease resistance, and production stability and quality (Meilleur and Hodgkin 2004, Engels 2006, Heywood 2007. They could contribute to a possible solution for the current scarcity and skyrocketing prices of natural vanilla. However, like many orchids, most of these wild vanilla species are under pressure, mainly due to deforestation, unsustainable collection from the wild and lack of conservation efforts (Verma et al. 2009, Wegier 2017. The genus Vanilla can therefore be considered as a group of flagship species for the effective protection of its CWRs, via an integrated approach of both conservation and sustainable cultivation Mosquera-Espinosa 2018, Herrera-Cabrera 2018). We believe that this can be achieved through an innovative concept of joint land sparing and land sharing (SPASHA; Fig. 1). Land sparing sets aside some land for conservation and some for intensive cultivation, while land sharing combines the cultivation of crops with the conservation of native vegetation elements in the same areas (e.g., agroforestry) (Fischer 2014, Phalan 2018. The SPASHA approach proposed here combines these two strategies in a single territory, by protecting natural vanilla populations and their pollinators inside forests, while using some of the neighboring land for vanilla cultivation in agroforestry systems where the natural interactions are kept as intact as possible.
In this study, we return to the regions of origin of aromatic vanilla species within the Vanilla genus, more specifically Costa Rica, and hypothesize that species distribution models (SDMs) can identify suitable areas for both the conservation and sustainable cultivation of vanilla CWRs, enabling us to implement the proposed SPASHA approach. SDMs, also known as bioclimatic envelope models, ecological niche models, and habitat suitability models, generate geographical maps of a species' environmental suitability. Thanks to the growing availability of freely accessible georeferenced species records (e.g., Global Biodiversity Information Facility, Botanical Information and Ecology Network) and environmental data (e.g., WorldClim, ISRIC), SDMs are now widely used in many ecological applications, such as the selection of priority areas for conservation, the prediction of climate change effects on species distribution ranges, and the determination of species' invasion risks (Guisan and Zimmermann 2000, Pearson and Dawson 2003, Guisan and Thuiller 2005, Phillips et al. 2006, Rodr ıguez 2007, Elith et al. 2010, Jim enez-Alfaro 2018. We focus on four wild aromatic Vanilla spp. (V. hartii, V. odorata, V. pompona, and V. trigonocarpa), of which aromatic profiles have been made for the species V. odorata and V. pompona (Soto 1999, Maruenda et al. 2013, and provide spatially explicit recommendations for the application of the SPASHA approach within our study area Area de Conservaci on Osa (ACOSA) in southwest Costa Rica, where the current main economic activities are ecotourism and agriculture focusing on cattle and palm oil (Sierra andRussman 2006, Hunt 2015). Using the national biological corridor network of Costa Rica (Sistema Nacional de Areas de Conservaci on 2018), we define corridors with high vanilla suitability that can be reforested into productive agroforestry systems with vanilla as (one of) the cash crop(s). This promising strategy can be integrated in forest landscape restoration programs where agroforestry systems provide an income for local communities, while also connecting the surrounding protected forest fragments through a corridor network harboring wild vanilla populations. Our study demonstrates that SDMs are a useful technique to implement theoretical concepts, such as the proposed SPASHA approach, in land management and policy. We believe that vanilla CWRs offer a high cultivation potential within sustainable agroforestry systems and could contribute to the improvement of current crop production systems, restore degraded landscapes, and enhance socio-economic activities while satisfying global food demands. The proposed model can be used in other areas where vanilla CWRs are present.

Study area
Costa Rica is located within the native growing ranges of several aromatic vanilla species and has a tradition of small-scale vanilla cultivation, focusing on the commercial species V. planifolia and more recently the hybrid of Costa Rica [V. planifolia 9 V. pompona] 9 V. planifolia (Quir os 2010). Costa Rica's lowland rainforests harbor several wild vanilla species of great interest for their potential in breeding programs, but understudied in terms of their biology and ecology (Azofeifa-Bolaños 2017). An important area regarding the presence of vanilla CWRs is the Area de Conservaci on Osa (ACOSA) in the southwest of Costa Rica (Fig. 2). ACOSA is an important endemism area and a biodiversity hotspot containing 2.5% of the world's marine and terrestrial biodiversity (Kohlmann 2010). Protected areas represent a third of the total land surface of ACOSA. The average annual precipitation ranges from 2500 to 6000 mm, with most precipitation occurring during the rainy season from May to November. The average annual temperature is 25°C, with local variations as a result of topography and other geographic differences (Holdridge 1967, Kappelle 2003. The dominant agricultural land uses in these areas are plantations of Elaeis guinensis (oil palm), Tectona grandis (teak), and Gmelina arborea (gmelina) and cattle grazing lands.

Study species
Vanilla Plum. ex Mill. is a pantropical genus belonging to the orchid family (Orchidaceae), comprising about 120 species of hemi-epiphytic and epiphytic climbing vines with branching stems Dressler 2010, Cameron 2011). Within the neotropical region approximately 60 vanilla species are now recognized, belonging to three distinct phylogenetic lineages: (1) the membranaceous clade that is sister to the other members of the genus, with 15 species (Vanilla subgen. Vanilla); (2) four Antillean species that are phylogenetically related to African species; and (3) the aromatic clade (Vanilla subgen. Xanata sect. Xanata) with 41 species. The commercial species V. planifolia, as well as several CWR species with aromatic fruits, fall within the latter clade (Soto Arenas and Dressler 2010). Little is known about the exact distribution of these species, and it is likely that nomenclatural changes are required in this taxonomically challenging genus (Soto-Arenas and Cribb 2010). Presently, 14 vanilla species have been reported for Costa Rica (V. costaricensis, V. dressleri, V. hartii, V. helleri, V. inodora, V. insignis, V. karen-christianae, V. odorata, V. phaeantha, V. planifolia, V. pompona, V. sarapiquensis, V. sotoarenasii, and V. trigonocarpa) in addition to the hybrid species mentioned above (Soto Arenas and Dressler 2010, Azofeifa-Bolaños 2017, Karremans and Lehmann 2018. Seven of these species (V. dressleri, V. hartii, V. karen-christianae, V. odorata, V. planifolia, V. pompona, and V. trigonocarpa) are found within our study region ACOSA (based on forest inventories carried out within ACOSA), all belonging to the aromatic clade.

Species distribution modeling
We used an ensemble modeling approach (i.e., a combination of different modeling algorithms) to model the distribution of vanilla crop wild relatives under current climate conditions and to identify suitable areas for vanilla conservation and cultivation. The ensemble modeling was implemented in the R package BiodiversityR (Kindt 2018), using the following twelve commonly used distribution modeling algorithms: Maxent (MAXENT), random forests (RF), support vector machines (SVM), flexible discriminant analysis (FDA), multivariate adaptive regression splines (MGCV), domain algorithm (DOMAIN), and a stepwise and non-stepwise implementation of boosted regression trees (BRT), generalized linear models (GLM), and ❖ www.esajournals.org generalized additive models (GAM). These algorithms are presence-absence algorithms, with the exception of Maxent (presence-background) and DOMAIN (presence-only). Both suitability and presence-absence maps were generated for Costa Rica and cropped to the extent of our study region ACOSA. This is further explained below in Distribution modeling.
Presence data and modeling extent.-A dataset of georeferenced presence locations of vanilla species occurring in Costa Rica was compiled from different sources, including records from online databases such as the Global Information Facility (Global Biodiversity Information Facility 2019), the Botanical Information and Ecology Network (BIEN; Enquist 2016) and Tropicos (tropicos.org 2019), scientific articles and reports (among others Soto Arenas and Dressler 2010, Azofeifa-Bolaños 2017, Karremans and Lehmann 2018), and forest inventories carried out within our study area and the rest of Costa Rica by ourselves, local universities (Universidad de Costa Rica and Universidad Nacional de Costa Rica), and partner organizations such as Osa Conservation and Ministerio de Ambiente y Energia-Sistema Nacional de Areas de Conservaci on (MINAE-SINAC). When using SDMs, the extent of the modeling range is of particular importance. It has been recommended to use occurrence data from the complete distribution range Area de Conservaci on OSA (ACOSA) in southwest Costa Rica. Corcovado and Piedras Blancas are the two terrestrial national parks, and the Golfo Dulce reserve is the forest reserve connecting these two national parks. Marino Bellena is the marine national park. Terraba-Sierpe is a protected wetland. Golfito refuge is a wildlife refuge bordering Piedras Blancas national park. The dotted lines represent the biological corridors as part of the National Bio-Corridor Program (PNCB) of Costa Rica, with the Osa corridor laying within our study region ACOSA. of the species. Two important reasons are as follows: (1) When taking a subset, it likely does not include the full environmental variation under which a species is occurring, and (2) many species are rare and are represented by just a few records (Reddy and D avalos 2003, Schulman et al. 2007, Beaumont et al. 2008, Hortal 2008, S anchez-Fern andez et al. 2011, Barbet-Massin 2012, Raes 2012, Gomes 2018. Completely opposite, others have suggested to select a restricted distribution range because of the occurrence of subspecies within for example geographically isolated areas and pooling these subspecies for analyses might mask biologically relevant spatial variation within the distribution range of a single species (Gonzalez et al. 2011, Gomes 2018. We chose to use the neotropics (Mexico toward south Brazil) as our extent for occurrence to prevent the above-mentioned deficiencies associated with a narrow geographical and thus environmental selection, and especially because of the scarcity of available georeferenced data on vanilla species (Pupulin and Karremans 2017). The compiled data were cleaned by removing all records with missing locality information, and all records with coordinates located in the ocean, coordinates with latitude or longitude equal to zero, or with latitude equal to longitude (Boyle and Smith 2010, Maldonado 2015, Zizka and Antonelli 2015. Four vanilla CWRs (V. hartii, V. odorata, V. pompona, and V. trigonocarpa) had more than 30 species occurrence records, which we considered enough to build accurate distribution models (Wisz 2008, Mateo et al. 2010, van Proosdij 2016, and all four of them are present within our study area ACOSA. The SDMs were carried out with minimum half of the presence and absence data within Costa Rica. If there were more data outside Costa Rica, a random subsample was taken outside Costa Rica with the same amount of data as inside. Background/pseudo-absence data and target group background approach. -In addition to the presence data, all modeling algorithms used in this study required the input of absence or pseudo-absence data except for MAXENT (presence/background) and DOMAIN (presence-only). Pseudo-absence and background points were selected using the target group approach (Phillips 2009), a method that has been used to improve the predictive performance of SDMs in the presence of spatially biased presence data (Phillips 2009, Elith et al. 2010, Gomes 2018. In the target group approach, pseudo-absence or background records are selected from grid cells with presence data of species that belong to a similar group of species as the target species, under the assumption that these locations reflect a similar bias as the sampling bias of the target species (Phillips 2009). In our case, the target group grid was constructed using the occurrence records of all hemi-epiphyte and liana species growing in Costa Rica (species list derived from tropicos.org). Presence data of these target group species were derived from GBIF and BIEN and cleaned using the same method as explained before. Background points were randomly selected from this grid. The pseudo-absence points were generated by selecting grid cells that did not contain records of the vanilla species being modeled and that had records of at least 20 other species. In this way, the pseudo-absence points were assumed to be relatively close to real absence points.
Predictor variables.-Data of 19 bioclimatic variables with a 1-km resolution were derived from the WorldClim database (Hijmans 2005). Furthermore, we selected eight continuous soil variables with a 250-m resolution from the ISRIC SoilGrid-s250 m database (Hengl et al. 2014), derived five topographic variables with a 250-m resolution from the Global Multi-resolution Terrain Elevation Data 2010 (GMTED2010; Danielson andGesch 2011, Amatulli 2018), and extracted cloud cover data with a 1-km resolution from the global remote sensing-derived cloud database (Wilson and Jetz 2016). We used the variance inflation factor (VIF) to avoid problems with multicollinearity of the predictor variables, using a stepwise approach that results in a set of variables that all have a VIF lower than 10 (Garc ıa 2015). Our final set contained nineteen predictor variables: mean diurnal range, maximum temperature of warmest month, mean temperature of coldest quarter, annual precipitation, precipitation seasonality, precipitation of warmest quarter, precipitation of coldest quarter, available water capacity, cation exchange capacity, clay content, coarse fragments, organic carbon content, pH index measured in water solution, silt content, aspect, slope, topographic position index, and cloud cover. Variables with a higher resolution of 250 m were resampled so all predictor variables had a spatial resolution of 250 m.
Distribution modeling.-The SDMs were evaluated using the area under the receiver operating characteristic curve (AUC) measure. The AUC measure is a threshold-independent statistic used in ecological modeling to assess the capacity of a model to discriminate between positive (presences) and negative (absences) instances (Phillips et al. 2006). The use of AUC has been criticized mainly because the increase of the geographical extent, and thereby the environmental extent, in which pseudo-absences are selected, causes an overestimation of AUC values (Lobo et al. 2010). To avoid this, we selected pseudo-absence and background points within a convex hull around the presence records, extended with a buffer of 10% of the longest distance between presence records. We crossvalidated the models with presence and absence data arranged in spatially independent blocks, using the blockCV package (Valavi 2019). Presence and absence data were divided into 100 km wide squared blocks arranged in 8 cross-validation folds. If the distribution range of a species was not large enough to use 8 blocks, a lower amount of folds was used. Using the ensemble.tune function of the BiodiversityR package, the weights of the different algorithms used in the ensemble model were optimized during each cross-validation run, and the average of these optimized weights was used to make the suitability predictions of the final ensemble model. Prior to making these predictions, the algorithms included in the ensemble model were calibrated again with all available occurrence data. To convert the continuous suitability predictions to presence-absence maps, the sensitivityspecificity criteria were used, one of the two methods recommended by Jim enez-Valverde and Lobo (2007), whereby the difference between sensitivity (true predicted presences) and specificity (true predicted absences) was calculated for the same 100 threshold values and the one which minimized that difference was selected. The resulting presence-absence maps were masked by the aforementioned extended convex hull to exclude areas without any occurrence records.
Variable importance and response curves.-We evaluated the importance of each predictor variable and visualized the response of the species to the three most important predictor variables using the evaluation strip method (Elith 2005).
The importance of the different predictor variables in the ensemble model predictions was calculated using a randomization procedure that measures the correlation between the original predicted suitability values and the predicted values where the variable under investigation is randomly permutated. If the contribution of a variable to the suitability predictions is high, the suitability predictions are more affected by the permutation of the variable in question, resulting in a lower correlation. Therefore, 1correlation was used as a measure of variable importance (Guisan et al. 2017).
Identifying vanilla conservation and cultivation areas: SPASHA approach.-Maps showing the national biological corridor network of Costa Rica (SINAC 2018) were overlaid with the generated presence-absence (1/0 values) maps to identify corridors with potential for the implementation of the SPASHA approach (i.e., areas with value 1 on the presence-absence maps). Focus for further specific land use recommendations was laid on the Osa corridor and its surroundings within our study region ACOSA. The presence-absence maps were cropped to the extent of ACOSA and overlaid with a recent land use map of this area (Shrestha et al. 2019, unpublished data) and a map depicting the protected areas (sinac.go.cr). Based on these final maps, we could make recommendations regarding land use management related to the SPASHA approach focusing on vanilla CWRs and based on the parameters given in Table 1. For example, if an area shows to be suitable for a certain vanilla species to occur based on the distribution modeling but is currently a forest patch laying outside a protected zone, it is seen as a priority area for vanilla conservation. Areas modeled to be unsuitable for a vanilla species to occur are left out from the Table 1, but are presented in the maps. All analyses were performed with R (version 3.4). We used the package BiodiversityR for ensemble suitability modeling and the packages ggmap, ggplot2, rgdal, rgeos, maptools, dplyr, tmap, mapview, leaflet, and raster for mapping operations.

Species distribution modeling
The ensemble species distribution models had a mean AUC value of 0.89. The mean AUC values for the four wild vanilla species (V. odorata, V. pompona, V. hartii, and V. trigonocarpa) were 0.86, 0.84, 0.94, and 0.93, respectively. Overall, the ensemble model, support vector machines (SVM), and maximum entropy model (MaxEnt) showed the highest performance among all species, as shown in Table 2. For all species, the ensemble model performed better than the individual algorithms. The most important predictor variables explaining the distribution of these four vanilla species are described in more detail in the paragraph below. The suitability maps representing areas with values ranging from 0 (unsuitable) to 1000 (highly suitable) and the AUC values of the individual modeling algorithms can be found in Appendix S1. Table 3 gives an overview of the contribution, with values ranging from 0 to 1, of each individual predictor variable included in the SDMs as well as the total contribution of soil, precipitation, temperature and topography separately, and the overall contribution of all predictor variables together. The most important predictor variables for V. odorata were the precipitation variables annual precipitation (0.31) and precipitation seasonality (0.24). The soil variables clay content (0.07) and coarse fragments (0.05) followed, but played a smaller role. The main variable explaining the distribution of V. pompona was the soil variable available water capacity (0.21), followed by maximum temperature of the warmest month (0.07) and mean temperature of the coldest quarter (0.05), but to a lower extent.
The key predictor variables for V. hartii were precipitation of the warmest quarter (0.29) and mean diurnal range (0.12). Additionally, soil variables such as coarse fragments (0.08), soil clay content (0.07), and available water capacity (0.06) were also fairly contributing to its distribution. The most important variables for V. trigonocarpa were maximum temperature of the warmest month (0.19) and annual precipitation (0.12). On average, precipitation and/or temperature variables (climate) contributed most to the distribution of the species V. odorata, V. hartii, and V. trigonocarpa whereas soil and also precipitation were more important for V. pompona. Overall, the predictor variables included in the model explained a considerable proportion of the vanilla species distribution. All the included predictor variables together contributed for 0.79, 0.55, 0.99, and 0.60 in the final distribution maps of V. odorata, V. pompona, V. hartii, and V. trigonocarpa, respectively.  In general, the mountainous areas (>1000 m) are less suitable or unsuitable for the studied vanilla species. The lowland areas (0-500 m) showed highest suitability, with some variations among the species, which depended on local or regional climate, soil, and topographic variables as mentioned above. For example, the drier areas in northwest Costa Rica are suitable for V. pompona, but not for the other species. Both Pacific and Caribbean coast lines have suitable environmental conditions for the four vanilla species to occur, but V. hartii and V. trignocarpa have a more restricted distribution compared to V. odorata and V. pompona, where V. pompona seems to be the most robust species with the widest potential distribution along Costa Rica.

Maps for priority conservation and sustainable cultivation of vanilla CWRs (SPASHA approach)
For each of the four vanilla species, the modeled presence-absence maps, derived from the suitability maps, were overlaid with the national biological corridor map of Costa Rica to identify the corridors with suitability for the implementation of the SPASHA approach. This is illustrated for the species V. odorata (Fig. 3), whereas the results for the other vanilla species can be found in Appendix S1. Overall, the biological corridor Osa and its surrounding areas within our study region ACOSA contain suitable areas for all four vanilla species to occur.
After cropping these maps to the extent of our study region ACOSA and overlaying them with a recent land use map and a map of the protected areas of ACOSA, a final map was created for each species. The final maps highlight areas for conservation and areas recommended for sustainable cultivation, focusing on the biological corridor Osa as priority area but also its surrounding areas. Such a final map is illustrated for V. odorata (Fig. 4), and land use policy recommendations can be made using this map together with Table 1. The maps for the other species can be found in Appendix S1. In general, green areas are forests suitable for vanilla to occur, giving these areas a recommended vanilla conservation status, but the priority further depends on the location (inside or outside protected areas such as national parks). For example, areas classified green located within the national parks are protected and thus of less concern, whereas those located outside the national parks should receive high priority regarding vanilla conservation, especially the Osa corridor as part of the PNCB. The brown and yellow areas are plantations and grazing lands, respectively, recommended for the implementation of sustainable agroforestry systems (AFS) with vanilla as a cash crop, especially in areas where natural vanilla populations potentially occur in the nearby forests. The white areas are unsuitable for vanilla to occur, so with no recommendations regarding vanilla conservation or cultivation.

Vanilla distribution maps
This is the first study that explored the use of species distribution modeling to identify priority areas for conservation and sustainable cultivation of vanilla crop wild relatives. The created distribution maps for Costa Rica had a mean AUC value of 0.86, 0.84, 0.94, and 0.93 for the species V. odorata, V. pompona, V. hartii, and V. trigonocarpa, respectively, which we considered highly satisfactory to use for further application of the SPASHA approach. We made land use recommendations for our study region ACOSA, using a recent land use map and protected areas map. The final maps (Fig. 4 for V. odorata; Appendix S1 for other species) can now be used to implement our recommendations within the corridor Osa and its surrounding areas.
The SDMs showed that most of the lowland (<500 m) rainforests along the Pacific and the Caribbean coast of Costa Rica possess suitable environmental conditions for all modeled vanilla species. The drier area in the northwest of Costa Rica seems to be only suitable for V. pompona, which might be the result of its thicker stem and leaves giving it a higher tolerance toward droughts. The central mountainous region is, as expected, not favorable for these vanilla species, probably due to unsuitable climate conditions at higher elevations (such as colder temperatures and a less pronounced dry season needed for flowering). The predictor variables included in the species distribution model contributed to a noteworthy percentage of the distribution of all vanilla species. Other factors, not included in the model, are most likely the presence of soil mycorrhiza and other microorganisms. It is well-known that the family Orchidaceae has an important relationship with root fungi. A single orchid species can be associated with several mycorrhiza fungi, with different functional consequences for the plant (Ruibal 2017). This also applies for vanilla species, where it has been shown that different fungi species (among others Ceratobasidium, Tulasnella) play an important role in seed germination and plant growth (Porras-Alfaro and Bayman 2008, Cameron 2011, Flanagan and Mosquera-Espinosa 2016, Gonz alez-Ch avez 2018). There are also still a lot of questions regarding the pollination and seed dispersal of wild vanilla species, which could be important factors in explaining the distribution patterns of vanilla CWRs.
The set of vanilla species initially considered in this study included all aromatic vanilla CWRs native to Costa Rica and contained presence data among their whole natural distribution range within the neotropics, but most of them could not be included in the final analysis because not enough observations were available to construct accurate species distribution models. We defined a minimum sample size of 30 observations (based on among others Stockwell and Peterson 2002), proven to be sufficient to perform reliable SDMs. Nevertheless, the lack of presence data of vanilla CWRs confirms the difficulty but also the importance of inventories and developing SDMs for these endangered species and stresses the need for more data collection as well as the creation of aromatic profiles of wild vanilla species, especially considering its economic value and current endangered status.

Joint land sharing/land sparing concept (SPASHA)
We provide spatially explicit recommendations for (1) the conservation of the four studied vanilla CWRs by protecting the forests that show high Fig. 4. Map highlighting areas for conservation and sustainable cultivation of V. odorata within the ACOSA region. The green areas are forests suitable for vanilla to occur, giving these areas a vanilla conservation status, but the priority further depends on the location (inside or outside protected areas such as national parks). The brown and yellow areas are plantations and grazing lands, respectively, recommended for the implementation of sustainable agroforestry systems (AFS) with vanilla as a cash crop, especially in areas where natural vanilla populations (potentially can) occur in the nearby forests. The white areas are unsuitable areas for vanilla to occur, so with no recommendations regarding vanilla conservation or cultivation. suitability for these species to naturally occur, and (2) the sustainable cultivation of these species within environmentally friendly farm systems around these forests, in areas suitable for the species in question but are currently dominated by monoculture cash crops and grasslands. We propose this as the joint land sparing / land sharing approach (SPASHA) and suggest it as a possible strategy to balance biodiversity conservation and agricultural production. Land sparing and land sharing have been projected as opposing strategies, where both social and biophysical factors determine which approach is most feasible or appropriate in a given landscape (Fischer 2011, Hulme 2011, Tscharntke 2012, Chandler 2013, Edwards 2014, Gilroy 2014. In the case of vanilla, we see a clear added value in the combination of both strategies, allowing the protection of natural forests where vanilla CWRs occur while surrounding, degraded lands can be used for the conversion into sustainable vanilla cultivation within environmentally friendly farm systems, such as agroforestry systems. For example, current grasslands or monocultures laying between forest fragments could be (partly) reforested for the creation of biological corridors where native trees species, who preferably provide economically interesting products, are planted and vanilla is cultivated as the cash crop. These agroforestry systems provide multiple ecosystem services and environmental benefits related to carbon sequestration, biodiversity conservation, soil enrichment, and air and water quality (Jose 2009). This approach allows local stakeholders in Costa Rica to get support from the Payments for Ecosystem Services Program (PESP), because (1) they reforest degraded areas for the cultivation of, in this case, vanilla within agroforestry systems and (2) conserve surrounding forests where natural vanilla populations occur (S anchez-Azofeifa 2007, Pagiola 2008, Farley and Costanza 2010. The latter is especially important because our SPASHA approach suggests, among other, the exchange of important natural interactions between forest and farmland. An example of such an interaction is animal pollination, a crucial activity for the reproduction of plants that depend on specific animals for their pollination, which is also the case for vanilla. Since the invention of the hand pollination technique by Edmond Albius in 1841, vanilla flowers in Madagascar's plantations have been pollinated manually and this successful technique is now applied worldwide (Havkin-Frenkel andBelanger 2010, Cameron 2011). So far, there is evidence from own field observations and a few previous studies (Lubinsky 2006, Pansarin and Pansarin 2014, Pansarin and Miranda 2016) that both the commercial species Vanilla planifolia and wild vanilla species within the subgenus Xanata, section Xanata, are naturally pollinated by orchid bees of the Euglossini tribe (Apinae). Very little is known about the pollinator species or their pollination success, but observations demonstrate that the bee species responsible for pollination vary upon vanilla species, probably due to differences in flower size, posing a reproductive barrier between species by selecting compatible pollinators. In general, it is known that bees provide important ecosystem services to agriculture, including pest control and crop pollination, and agricultural landscapes can in turn provide habitats for several bee species (Bianchi et al. 2006, Klein 2009, Kennedy 2013, Martins et al. 2015. This important group of animals is threatened by the loss of natural and semi-natural habitats, extensive monoculture plantations, and increased pesticide and herbicide use (Freitas 2009). In the case of vanilla, we think that our proposed SPA-SHA concept can protect both CWRs and CWRpollinator as well as other interspecific interactions. However, further study on natural pollination of both commercial vanilla species and vanilla CWRs, as well as flower traits and aromatic flower compounds in relation to pollinator species, is highly needed. In this way, intra-and interspecies compatibility can be determined and affect the success of the SPASHA approach.
Crop wild relatives have long been recognized as priority species for conservation because their extinction threatens the genetic base of many of the world's crops (Heywood 2007, Castañeda-A lvarez 2016. Their protection started with ex situ conservation programs, which were not as successful as expected, especially for species that are characterized by small and dispersed, genetically distinct populations with low seed viability and complex interspecific relationships within their natural habitat, which is the case for vanilla species (Alomia 2017, Azofeifa-Bolaños et al. 2018. High maintenance costs of clonal collections and problems regarding the regeneration of stored material and seed recalcitrance pose additional challenges for the ex situ conservation approach (Pritchard 2012, Braverman 2014. Therefore, interest in in situ conservation (Harlan 1992) has increased over the years and is now considered a key element to conserve CWRs. Joint land sparing and land sharing are an example of in situ conservation (the conservation of forests with vanilla populations) and ex situ conservation (the development of areas where CWRs will be cultivated). The starting point of CWR conservation are largescale inventories (Gadgil 1993). In 1985, Bioversity International established a list of genera for priority surveys and in situ conservation (IBPGR), stimulating country's interest and action toward CWR inventories. Hodgkin et al. (2001) stressed the priority of endangered species given the potential negative impacts of climate change and habitat fragmentation on unstable and isolated populations (Wilcove and Chen 1998, Brigham et al. 2002, Pitman 2002. Many species within the orchid family are currently facing extinction risks. Despite substantial conservation emphasis on rare orchids, populations continue to decline (Gale 2018). Almost all the species within the genus Vanilla that have been assessed in the IUCN Red List are listed as endangered , Wegier 2017. Their natural populations are threatened by deforestation and land degradation as well as uncontrolled harvest from wild populations (Schl€ uter et al. 2008, Flanagan andMosquera-Espinosa 2016). Since vanilla is a high-value crop facing many problems in the current vanilla production systems and given that its CWRs hold beneficial features for crop improvement, we suggest the listing of its CWRs and the start of large-scale inventories along its distribution range as high priority actions. Especially because these results can be included into management plans and implemented on broader scales.

Implementation in land use policy
The study emphasizes the proposed management interventions on the biological corridors as part of the National Bio-Corridor Program (PNCB) and surrounding forest fragments. The PNCB, a network of connecting structures between state conservation areas, was created in 2006 and reformed in 2017 as a strategy of biodiversity conservation. The goal of this program is to promote the conservation and sustainable use of the biodiversity in Costa Rica, from the perspective of a functional and structural ecosystem connectivity (SINAC 2018). Up to present, a strategic plan of the protection and management of these corridors in collaboration with local stakeholders (Plan Estrat egico 2018-2025, as part of PNCB) is formed but no action has been taken yet. There is a plan for the establishment of incentive systems and alternative income-generating opportunities within these corridors, where the involvement of local parties is key (SINAC 2018). As one of the first in contributing to this plan, we propose the implementation of the SPASHA approach within these corridors and their surroundings through the creation of diverse and productive landscapes with vanilla CWRs production as a high-value certified crop on degraded lands while protecting the surrounding natural forests. Besides the corridor Osa, we also consider the suitable areas within the Golfo Dulce Forest Reserve and its surroundings as areas of importance. Although defined as a forest reserve, personal annotations and observations from MINAE-SINAC indicate that this area is still exposed to deforestation practices for land conversions and illegal logging. With our approach, we could provide sustainable land use practices within these areas to prevent this from happening. It contributes to one of the goals of the PNCB, to give these areas an environmental as well as economic value, and offering alternative income generation for the local communities.

CONCLUSIONS AND FURTHER RESEARCH
Costa Rica is recognized for its positive attitude toward biodiversity conservation strategies and natural resource management in general, with clean energy production, the creation of 11 conservation areas, an extensive national park system, and a national biological corridor network as prime examples. However, unprotected lands are still being converted into rather unsustainable land use practices, with the main agricultural practices focusing on monocultures of oil palm, pineapple, banana, and teak or grasslands for cattle grazing. Moreover, it is the leading country regarding the use of biocides per unit area (Faostat 2019). Therefore, Costa Rica serves as an interesting country to apply the joint land sparing/ land sharing (SPASHA) strategy. Given the growing awareness of the economic potential of vanilla species, they should receive priority action for conservation and sustainable use, recently pointed out by among others Flanagan and Mosquera-Espinosa (2018) in Colombia and Herrera-Cabrera (2018) in Mexico. We propose the incorporation of the SPASHA strategy within land management plans in order to reduce the increasing pressure of unsustainable extraction from natural populations and to promote sustainable cultivation practices. For this study, we focused on the corridor Osa because of the availability of a recent land use map of our study region ACOSA. Future research should involve the creation of recent, high-quality land use maps of the whole country, for the implementation of our findings regarding land management plans within the biological corridors and their surroundings at a national level. We also suggest to continue collecting records on vanilla CWRs through forest inventories and collaboration between countries within the native ranges of aromatic vanilla species. In this way, the SPASHA approach could be implemented on a broader scale and we can collect more data on distribution, taxonomy, biology, ecology and potential use of different vanilla CWRs along the neotropics. By using modern technologies such as iNaturalist, local communities can get involved in the gathering of georeferenced data as well as in the collection of data on flowering, fruiting, and specific habitat preferences. Furthermore, we did not model the potential niche shifts of vanilla species under climate change but suggest this should be explored in a further study. The latter should include data on the distribution of associated pollinators, seed dispersers, microbial communities of leaf litter and soil, and mycorrhizae.