Climate change impact assessment on the potential rubber cultivating area in the Greater Mekong Subregion

In order to map potential shifts of rubber (Hevea brasiliensis) cultivation as a consequence of the ongoing climate change in the Greater Mekong Subregion (GMS), we applied rule-based classifications to a selection of nine gridded climatic data projections (precipitation and temperature, and global circulation models (GCMs)). These projections were used to form an ensemble model set covering the representative concentration pathways (RCPs) 4.5 and 8.5 of the Fifth Assessment Report of the Intergovernmental Panel on Climate Change at three future time sections: 2030, 2050 and 2070. We used a post classification ensemble formation technique based on a majority outcome of the classification to not only provide an ensemble projection but also to spatially track and weight the disagreements between the classified GCMs. A similar approach was used to form an ensemble model aggregating the involved climatic factors. The level of agreement between the ensemble projections and GCM products was assessed for each climatic factor separately, and also at the aggregate level. Shifting zones with high confidence were clustered based on their land use composition, physiographic attributes and proximity. Following the same ensemble formation technique and by setting a 28 °C threshold for annual mean temperature, we mapped areas prone to exposure to potentially excessive heat levels. Almost the entire shift projected with high certainty was in the form of expansion, associated with temperature components of climate and temporally limited to the 2030 time window where the total area conducive to rubber cultivation in the GMS is projected to exceed 50% by 2030 (from 44.3% at the turn of the century). The largest detected cluster (41% of the total shifting area), which also is the most ecologically degraded, corresponds to Northern Vietnam and Guangxi Autonomous Region of China. The area exposed to potentially excessive heat is projected to undergo a 25-fold increase under RCP4.5 by 2030 from 14568 km2 at the baseline.


Introduction
Natural rubber is a key industrial commodity with wide application in manufacturing of a very diverse range of products. Although rubber-bearing plant species such as Taraxacum kok-saghyz and Parthenium argentatum have lately reemerged on the research and development scene as potential alternative sources of natural rubber (van Beilen and Poirier 2007a, 2007b, Rasutis et al 2015, Kreuzberger et al 2016, Dong et al 2017, Ramirez-Cadavid et al 2017, Soratana et al 2017, the Para rubber tree (Hevea brasiliensis) has retained its status as the sole viable source of natural rubber, which does not seem set to change in the near future (Cornish 2017). Global consumption of natural rubber has exceeded 12 million metric tons in the last three years according to the International Rubber Study Group (IRSG 2017). Raising demand has been matched and to some extent surpassed by increases in production. Global trends of natural rubber production and consumption, and the harvested area are illustrated in figure 1.
Recent decades have been associated with expansion (and to some extent shift) of rubber cultivation Global trends for consumption, production and area under rubber cultivation. Segmented regression lines reveal the shifts in trends: 1996 is the estimated year before which a 118.9 thousand ton increase per year explained the growing consumption trend, accelerating to 220.4 thereafter, while for production, the slope has shifted from 122.3 to 304.9 thousand tons per year by 1998 and year 2002 appears to be the most efficient breakpoint explaining the increasing trend of the global area under rubber cultivation, surging from 89.3 to 287.9 thousand hectares added each year. We have used R package 'segmented' (Muggeo 2003) version 0.5-1.4 to generate this figure from FAOSTAT (production and area) and IRSG (consumption) data (FAOSTAT 2017, IRSG 2017. Inkscape 0.91 was used for visual optimization. zones from the traditional rubber growing regions (the 10 • S to 10 • N equatorial belt) to higher latitudes and longitudes (Priyadarshan et al 2005, Ziegler et al 2009, Li and Fox 2012, Ahrends et al 2015, Chen et al 2016a, Chen et al 2016b. Thailand, the leading rubber producing country since 1990, which also has had the highest share of the global area converted each year to rubber cultivation (30% on average since the turn of the century), can well illustrate the situation (figure 2).
The Greater Mekong Subregion (GMS) is an economic cooperation program consisting of six nations: China (Yunnan province and Guangxi autonomous region), Vietnam, Thailand, Laos, Myanmar and Cambodia. The GMS covers more than 2.5 million km 2 of mainland Southeast Asia (MSEA), about 84% of which overlaps with the Indo-Burmese mega biodiversity hotspot (Myers et al 2000, Mittermeier et al 2004. It stands for a substantial share of global rubber production (46.7% in 2014) 3 almost exclusively coming from monocultures. Since its inception in the early 1990s, the GMS has in general, and its formerly isolated members (Myanmar, Laos and Cambodia) 3 This figure is mainly based on FAOSTAT data. As two Chinese provinces of Hainan and Guangdong contribute to the Chinese national production, their share (46.2% in 2014 as mentioned in the China Statistical Yearbook 2016 www.stats.gov.cn) has been deducted. In case of Laos for which FAO data is not available, United Nations Commodity Trade Statistics Data (comtrade.un.org) was used in combination with the historical commodity prices (www.indexmundi.com) to estimate the national rubber production: 56 thousand tons. in particular, been undergoing rapid socio-economic change through regional development programs and transboundary investments in all conceivable sectors. At the same time, ecological degradation through accelerated landscape transformation has been observed. Heavy expansion of rubber monocultures and their spread to new areas have had a notable contribution to deforestation, habitat fragmentation and biodiversity loss (Li et al 2007, Ahrends et al 2015, He and Martin 2015, Häuser et al 2015.
In response to concerns about the ecological implications of the rapid expansion of rubber monocultures mostly replacing forests and swidden agriculture in MSEA, remote sensing techniques are regularly used to monitor land use conversion to rubber cultivation (e.g. Li and Fox 2011a, 2011b, 2012, Dong et al 2012, Senf et al 2013, Fan et al 2015, Grogan et al 2015, Chen et al 2016a, Kou et al 2017. More recently, remote sensing has been used to track additional details such as the rubber plantation age (Koedsin and Huete 2015, Kou et al 2015, Beckschäfer 2017, Trisasongko 2017. Climate is one of the defining factors of the potential geographic extent for cultivation of any crop, and Para rubber is no exception. Momentous ongoing change in Earth's climate attributed to human activity (Collins et al 2013, Kirtman et al 2013, Lewandowsky et al 2016, Thorne 2017, Medhaug et al 2017, Berger et al 2017 is comprehensively acknowledged by the scientific community (Cook et al 2016). Some forecasts of the future potential geographical range for Para rubber in differ-  (NSO 1994) and the agricultural statistics yearbooks of Thailand (2009 and data (available at www.oae.go.th) and the GADM administrative division shapefiles (2.8) were used. Maps were generated in ArcGIS 10.2.2 and visually optimized in Inkscape 0.92. ent parts of MSEA, mainly based on ecological niche modeling (Ray et al 2014, 2016a, Ahrends et al 2015, Liu et al 2015 and bioclimatic stratification (Zomer et al 2014) have recently been published.
Gridded data of climatic factors simulating likely future conditions are essential inputs for forecasts. Global circulation models (GCMs) are useful sources of information commonly exploited to assess the potential impacts of climate change. Various institutions are engaged in creating such datasets and provide dozens of potential choices as input. Variations among GCMs, which mainly rise from structural and parameterization differences (Semenov and Stratonovitch 2010), can help to provide a means to capture and explore some of the projection uncertainties that have to be accounted for in order to obtain a realistic and scientifically sound image. Variabilities observed in sets of comparable simulations prompt some key choice questions, starting with whether a single simulation would suffice or a multi-member ensemble is needed for a reasonably robust forecast. In the latter case, can using the largest possible ensemble be a legitimate decision or could a reduced set of simulations perform better while minimizing the computational cost? Based on what criteria should a shortlisting take place? Should an average of all set members be used as ensemble or (considering the spatial nature of the data) is there a better option? How should the uncertainties (dispersion) inherent in the input differences (an important but so far overlooked factor) be handled? And how should these uncertainties be communicated in a comprehensive and useful way?
Potential phytosanitary deficiencies as well as growth and yield failures due to crop exposure to excessive levels of ambient temperature are some of the more unsettling aspects of climate change. Despite the existing evidence for this matter (Abd Karim 2008, Kositsup et al 2009, Yu et al 2014, Golbon et al 2015, Jayasooryan et al 2015, Nguyen and Dang 2016, setting a clear-cut threshold for heat stress is still a debatable subject.
Here, we apply rule-based geographical classification to a selection of the downscaled IPCC AR5 climatic projections in order to map the potential geographical zones projected to be climatically suitable for Para rubber cultivation, or exposed to excessive heat, in MSEA in three time sections centered on 2030, 2050 and 2070 while accounting for and presenting the classification uncertainty.

Data
We used the WorldClim dataset (version 1.4, Hijmans et al 2005) to generate the baseline climatic map and  (2015), which ranked IPCC AR5 GCMs according to their regional performances and recommended a subset of eight to ten GCMs, avoiding the least realistic models while retaining the maximum plausible dispersion. Nine GCMs were selected using the regional plausibility rankings: The GCM data were provided by the Climate Change and Food Security Program data portal of the International Center for Tropical Agriculture (available at www.ccafs-climate.org) and were downscaled to 30 arc sec (∼1 km) resolution using delta method (Ramirez-Villegas and Jarvis 2010). Two of the four main climate change scenarios recognized by the IPCC AR5 were considered in this study: Representative Concentration Pathways (RCPs) 8.5 and 4.5. RCP 8.5 is a high greenhouse gas (GHG) emission scenario comprising no stabilization of the atmospheric GHG concentrations leading to 8.5 W m −2 of radiative forcing by 2100 and a globe over 4 • C warmer than the pre-industrial era. RCP 4.5 is a moderate scenario accommodating GHG concentration stabilization by 2070 and radiative forcing of 4.5 W m −2 (2.5 • C temperature rise) by the end of the 21st century (Riahi et al 2011, Thomson et al 2011. Land use data (see supplementary material figure S1 available at stacks.iop.org/ERL/13/084002/mmedia, Hoskins et al 2016), the biodiversity intactness index (BII) created by Newbold et al (2016, supplementary material figure S2) and the USGS GTOPO30 digital elevation model were used to cluster and describe the potential future expansion/retraction zones. We have also used the administrative divisions (GADM) shapefiles (available at www.gadm.org) in this study. Steps involved in intra-annual temperature distribution suitability classification (as an illustration case). Continuous monthly mean temperature gridded data (RCP 4.5 for the 20 year period centered on 2030) were used to generate binary layers by setting a stress threshold (23 • C). All 12 binary layers originating from the same GCM were summed to produce the layers reflecting the number of months projected to be below the threshold (abbreviations AC to NO denote the corresponding GCMs). These layers were reclassified to three levels: optimal, suboptimal and prohibitive . The ensemble classification map was generated by extracting the majority outcome of all GCMs for each grid cell . The uncertainty layer reflects the consensus level among GCMs leading to the ensemble and was produced by counting the number of GCMs participating in the formation of the majority for each given grid cell . Panel shows the geographic extents of the frame selected for illustration. All layers used in each step were assigned equal weights and the arrow color difference is only for visual clarity. ArcGIS 10.2.2 and Inkscape 0.91 were used for the generation of this figure.

Methods
Five climatic suitability criteria adapted from Rivano et al (2015) listed in table 1 were used in this study. As mentioned by Thompson et al (2013) and Stephens et al (2012) it is essential to avoid averaging for ensemble formation as it leads to information loss on variation.
Here we conducted the complete classification process on the involved gridded variables separately for each GCM and formed the ensemble product by the majority outcome for each grid cell overlaid with a simple uncertainty measure reflecting the strength of the majority. The total annual precipitation and the mean annual temperature layers were directly categorized to optimal, suboptimal and prohibitive ranges for each GCM, time section and scenario. The ensemble suitability projections were generated for each 'criterion × time section × scenario' combination consisting of the suitability class returned by the majority of the GCMs for every grid cell and a corresponding uncertainty layer reflecting the strength of the consensus on the class assigned to each ensemble grid cell ranging from full agreement (9/9) to mere majority (5/9). The monthly mean temperature and the monthly precipitation gridded data went through a similar process with two additional steps (see figure 4), summarizing the intra-annual distribution of precipitation and temperature.
By overlaying the classification outcomes of the climatic layers, each grid cell was assigned one of the following summarizing classes: 'AllOpt' where all climatic layers returned an optimal classification, 'SubOpt' where at least one layer was described as suboptimal and none as prohibitive, 'SingProh' where only one layer was in prohibitive range and 'MultProh' with more than one climatic criterion in the prohibitive range. The aggregate uncertainty layers were also overlaid to produce an aggregate uncertainty layer in four levels: (1) full agreement among GCMs for all four criteria, (2) only one criterion projected with 7 or 8 from 9 majority (and all other criteria possessing stronger consensus), (3) only one criterion projected with 5 or 6 from 9 majority (and stronger consensus in all other criteria) and (4) two or more criterion projected with 5 or 6 from 9 majority.   Rivano et al (2015). The number of months with mean temperature below 23 • C is referred to as intra-annual temperature distribution and the number of months with precipitation below 50 mm as intra-annual precipitation distribution.
A point shapefile representing grid cells in the raster data was created for the shift zones with high aggregate certainty (levels 1 and 2) to which the corresponding land use, BII, altitude, slope, longitude and latitude values both in original and standardized form were extracted. We used the Grouping Analysis tool of the ArcGIS 10.2.2 to form clusters based on the standardized attributes and illustrated the outcome using 'ggplot2' 2.1.0 package in R.
Sankey diagrams are illustration tools suitable for description of multidimensional and hierarchical categorical data and are most often used to show material or energy flows through network systems. Geographical classification dynamics over time can also be very efficiently presented by Sankey diagrams. As demonstrated by Cuba (2015), Sankey diagrams are superior to crosstabulation matrices in reflecting land use dynamics, particularly when multiple time sections are of interest. We generated Sankey diagrams to illustrate the climatic suitability class shifts projected to occur under RCP 4.5 and RCP 8.5 for each adjacent pair of time sections using the D3.js JavaScript library developed by Bostock et al (2011).
Using the ensemble formation technique, we created an 'excessive heat' layer distinguishing the area associated with annual mean temperature exceeding 28 • C at the baseline and traced its potential expansion under the two RCPs overlaid with corresponding uncertainty layers. This criterion, however, was not used as an upper limit for transition to suboptimal or prohibitory conditions in the former steps.

Single criterion classification
Climatic conditions in the study area at the baseline and the ensemble projections for the four climatic criteria (separately classified) are presented in figures 5 and 6. Largest projected shifts (expressed as proportion of the total studied area) are observed for the annual mean temperature and the intra-annual temperature distribution moving from baseline to the 20 year time window centered on 2030. Considering the annual mean temperature, 25.13% of the total area (642416 km 2 from 2556370 km 2 ) is projected (21.79% projected with full GCM consensus) to migrate from prohibitive and suboptimal range to classes more conducive to rubber cultivation under RCP 4.5. The RCP 8.5 ensemble projection suggests this figure to be 28.55% (23.38% with full agreement). For intra-annual temperature distribution, 20.18% (16.59% with full agreement) of the total area is observed to experience such a transition under RCP 4.5 and 23.96% (17.67%) under RCP 8.5 from baseline to the 2030 period. Moving to 2050 and 2070 time periods, the emerging more suitable areas regarding the two aforementioned factors are of much smaller size and paired with higher degrees of uncertainty. The persistence of the new conditions in an area that has gone through climatic shift is relevant but not necessarily traceable in Sankey diagrams (figure 5). Considering annual mean temperature under RCP 4.5, 14.17% of the total area is projected with high certainty to remain in the new class after shifting from prohibitive to suboptimal or suboptimal to optimal range and 16.58% under RCP 8.5. For the intra-annual temperature distribution, these figures are projected to be 14.17% and 16.55% respectively. Unlike temperature components of climate, the projected shifts observed in precipitation components were bilateral, associated with low certainty (i.e. high disagreement among GCMs) and smaller in size. The largest area projected to experience shifts in the precipitation class by 2030 was observed for intra-annual precipitation distribution summing to 7.02% of the total investigated area (6.01% moving from prohibitive to suboptimal or suboptimal to optimal and 1.01% vice versa). 80.63% of it (equal to 5.66% of the total area) has been projected by mere majority (i.e. the lowest possible certainty level).
Comparing single criterion GCM projections with their corresponding ensembles (table 2) reveals that ACCESS1.0 has returned the closest single criterion maps to the ensemble with an average overlap of 92.42% across all 24 possible criterion-RCP-time period combinations followed by BCC_CSM1.1, MPI-ESM-LR and CCSM4.

Aggregate classification
The geographical and temporal dynamics of the projected climatic suitability classes at the aggregate level are illustrated in figures 7 and 8. The area projected to retain its aggregate climatic class across the investigated time span (by 2070) is projected to be 72.83% of the total area under RCP 4.5 and 66.23% under RCP 8.5. By the time window centered on 2050 these projections sum to 74.98% and 72.89% and by 2030 to 77.63% and 78.22% of the total area respectively. From the total projected class-shifting area by 2030, 26.78% (6.01% of the studied area) was projected with maximum certainty (i.e. full agreement among GCMs in all four criteria) under RCP 4.5 and 26.50% (5.77%) under RCP 8.5. It was projected to decline to 14.09% (3.53%) and 17.39% (4.71%) for the baseline to 2050 time period and further reduction to 9.49% (2.58%) and 7.10% (2.40%) for 2070 respectively. Performance similarity of single GCM aggregate classification maps with the ensemble is presented in table 3 where ACCESS1.0 returned the closest results to the ensemble.
Restricting the investigated time window to the 20 year period centered on 2030 and the area to where climatic conditions are projected with high (the upper two levels) certainty to shift from prohibitive to rubber cultivation to suboptimal or optimal, our projections detected 195928 km 2 (7.70% of the total investigated area) under RCP 4.5 and 238734 km 2 (9.38%) under RCP 8.5, which are presented in figure 8. Using grouping analysis we detected eight major clusters based on land use composition, physiographic attributes and proximity. Northernmost potential expansion was projected to verge on 27 • N of the Irrawaddy basin and the altitudinal limit to exceed 1400 m a.s.l. in clus- ). The overall baseline state of biodiversity in these clusters is presented in figure 10 using the BII. Steffen et al (2015) proposed a safe limit value of 0.9 (maximum 10% decline) for the BII. The largest cluster (cluster 1) corresponding to Guangxi Autonomous Region of China and Northern Vietnam is the most ecologically degraded and accommodates 92.47% of the area already below the safe threshold.

Exposure to excessive heat
Projected exposure to annual mean temperature levels exceeding 28 • C in the study area is presented in figure  11. Based on WorldClim data, the total baseline area with this characteristic is limited to 14570 km 2 (less than 0.6% of the total investigated area) located between 12.33 • N, 100 • E and 16.33 • N, 101.50 • E in Thailand. Ensemble projections based on 7/9 to 9/9 majority classification suggest that by 2030, under RCP 4.5 this area may increase 25 fold (14.3% of GMS) and 35 fold (20.5% of GMS) under RCP 8.5 stretching northwards to 22 • N in the central parts of the Irrawaddy basin. By 2050 however, this criterion may be associated with 23.2% of the total area under RCP 4.5 and 31.2% under RCP 8.5 increasing respectively by 2070% to 26.5% and 38.9%. Zomer et al (2014) conducted a study focusing on the potential changes in the area conducive to rubber cultivation in Xishuangbanna, Yunnan, China using environmental stratification while averaging all four AR5 RCPs, which suggested an increase from 33.5% to 74.5% of the total prefecture area by 2050. Our findings for the same temporal and spatial frame are 52.5% (43.7% with high certainty) under RCP 4.5 and 83.1% (60.1%) under RCP 8.5. Ray et al (2016bRay et al ( , 2016b used MaxEnt ecological niche modeling tool, exploring the rubber producing Western Ghats and the North-East regions of India, and noted a substantial attachment of the projection outcome to the region used for calibration. If Amazonia was used for model calibration, only a very limited southern part of Western Ghats was returned as suitable by MaxEnt while established rubber growing regions were left out. They observed the same limited transferability pattern while calibrating MaxEnt with each of two Indian rubber producing Figure 9. Areas projected with high certainty to become climatically suitable for rubber cultivation by 2030. Land use composition (Hoskins et al 2016), physiographic composition (USGS GTOPO30) and the BII (Newbold et al 2016) were used to group the parts of the study area that were projected with high ensemble consensus to become climatically suitable for rubber cultivation into eight clusters using the grouping analysis tool of ArcGIS 10.2.2. Violin plots (bottom panel) were generated using the 'ggplot2' 2.1.0 R package and visually optimized in Inkscape 0.92. The terms primary and secondary habitat represent 'undisturbed natural' and 'recovering, previously disturbed natural' habitats respectively. The variables shown above are adjusted to share zero mean and unit variance. For the original scale, please see supplementary material figures S4-S6.

Contrasts and conjunctions with comparable studies
regions projecting for the other, one at a time. They reached plausible projections only by pooled occurrence points for parameter estimation. Ahrends et al (2015) investigated the expansion trends of rubber cultivation in roughly the same geographical frame as this study and concluded that this land use is stretching into increasingly less suitable zones jeopardizing biodiversity and landscape functions. They included a typhoon damage risk assessment based on historical tropical cyclone tracks which, when compared with the area projected with high certainty in this study to become climatically conducive to rubber cultivation by 2030, suggests current typhoon risk zones to overlap only with parts of clusters one (13.2%) and three (2.2%). This overlapping area in cluster one is limited to a 50 km inland buffer of the Guangxi coastline between 106.50 • E and 109.66 • E. Recent studies on the influence of climate change on western North Pacific tropical cyclone tracks, however, project reductions in both frequency and intensity of typhoons in future for our area of interest, mainly due to northward diversion (Colbert et al 2015, Kossin et al 2016, Zhang and Wang 2017. Liu et al (2015) projected the change in the area with potential for Para rubber cultivation in China covering all five provinces with a rubber cultivation background (Hainan, Yunnan, Guangdong, Guangxi and Fujian) using ecological niche modeling and reported a 15% increase by 2050 from about 400000 km 2 in 2010.
With the exception of cluster 1, which encompasses a major biotically compromised (BII < 0.9) area share, most of the regions projected to gain climatic potential for rubber cultivation are chiefly composed of intact primary habitats (figure 10 and supplementary material figure S4). In these areas land use modifications of significant scale require serious attention to the potential impacts on the ecological integrity and ecosystem functions and services. The ongoing improvements in the scientific understanding and practice of concepts such as rubber based agroforestry systems (Langenberger et al 2017) and Green rubber eco-certification (Kennedy et al 2017) offer promising options for environmentally friendly rubber cultivation, particularly as support from the smallholder side for participation in ecosystem protection appears to grow (Min et al 2018, Wigboldus et al 2017.

Strengths and limitations of the projection approach
Hevea brasiliensis is not only a plant and therefore a sessile species, but also a crop subject to non-natural sources of influence (e.g. breeding and crop management), which may affect the reliability of species distribution models if based on biased presence and pseudo-absence records. From our point of view rule-based models tend to be less prone to circular reasoning but risk engaging non-accurate classification rising from mis-estimated or dated tolerance thresholds (e.g. due to breeding).
We chose to assign equal weights to the climatic criteria involved in this study, and also to the GCMs forming the ensemble single criteria layers. However, we acknowledge that a non-equal weight approach based on justified quantification of the influence associated with each criterion or its ensemble projection homogeneity (in the case of GCMs, based on data quality) is plausible.
Non-climatic factors (e.g. soil conditions, land physiography, labor and market access), which are known to be decisive in suitability for rubber cultivation, were not involved in this study. Coverage of a broad range of suitability determining factors in a single study faces serious technical challenges. Different variables can often not be processed with a general approach as the scale relevant for some factors may not necessarily match the scale suitable for the others. The availability and quality of data in a standardized form are also two crucial limiting features. However, some factors relevant in smaller scale (e.g. soil properties) can be nested in those relevant in larger scale (e.g. climatic conditions) by subsequent localized assessments. This requires the provision of the outputs of studies such as this in a modular form exploitable for third parties. The KMZ files accompanying this manuscript do not only provide the findings unchained from resolution loss, but can also be used by future studies as a base to expand upon.
Although recent trajectories of GHG emissions are closest to the RCP 8.5 (supplementary material figure S7), this climate change scenario incorporates some assumptions concerning the use of fossil energy resources that are in the long-run technically improbable (Capellán-Pérez et al 2016, Ritchie and Dowlatabadi 2017a, 2017b. In view of the concerns and evidence regarding rapid changes in land use and climate, it is counterintuitive to use the early years of the last decade as baseline. Nevertheless, most required underlying data components are being revised not with emphasis on updating but on resolution (e.g. Newbold et al 2016) or precision (e.g. Fick and Hijmans 2017).
Compared with the lower temperature tolerance limits known for Para rubber, upper thresholds and consequences of exposure to high levels of ambient temperature are not well understood. The global area already exposed to annual mean temperature above 28 • C (supplementary material figure S9) does not match typical rubber growing regions. In the case of the GMS, a comparison between figures 2 and 11 underlines this point. Mesocosm experiments (Stewart et al 2013, Bestion et al 2015, Fordham 2015 and other manipulation methods that have recently gained prominence in studies aiming at a better understanding of the responses of the organisms to a warming climate can illuminate the way for H. brasiliensis as well.
The methods developed in this article are applied to a relatively restricted case study, rubber cultivation in the GMS. Nevertheless, the potential for transferability to other world regions and other cropping systems is very high, as the vast majority of datasets used is freely available for scientific purposes. The phenological and physiological crop specific background data for other crop plants can be collected from text books and literature reviews. Potential applications that come to mind might be the potential suitability for palm oil plantation systems, coffee agroforestry or bio-economically important crops such as sugar cane and maize and their potential northern distribution limits.
In order to broaden the audience of this study and to facilitate the use of its outputs for potential decision makers, we have produced two KMZ files Table 2. Classification agreement between the single criterion climatic data simulations and their ensemble.
Resemblance of the climatic classification by each of the nine simulations used in this study with their ensemble is expressed as proportion (%) of the sum of the areas with matching classification to the total area. Color-code reflects five levels: below 75% , 75%-90% , 90.1%-95% , 95.1%-99% and above 99% .  Table 3. Classification agreement between the data simulations and their ensemble at the aggregate level.
Resemblance of the aggregate climatic classification by each of the nine simulations used in this study with their ensemble is expressed as proportion (%) of the sum of the areas with matching classification to the total area. Color-code reflects five levels: below 75% , 75%-80% , 80.1%-85% , 85.1%-90% and above 90% . Nine IPCC AR5 simulations of representative concentration pathways RCP 4.5 and RCP 8. (one for each RCP), which summarize the information behind figures 5, 7 and 11, covering the baseline and the 2030s time windows. These files can easily be loaded in Google Earth to check the conditions for a given location by clicking. The KMZ files and all of the figures in high resolution are available on zenodo.org (https://doi.org/10.5281/zenodo. 1312769).

Conclusion
Even though the climatic change in the GMS is projected to be predominantly in the direction of higher suitability for rubber cultivation, the expansion of the climatically optimal area is projected to be minimal. When including the exposure to annual mean temperatures exceeding 28 • C (current estimate of excessive heat for Hevea rubber) as a limiting factor, then even a heavy reduction in the total climatically optimal area is likely to occur (see figure 8).
Across the time span investigated in this study (limited to 2070), about half of the new area with climatic potential for rubber cultivation is projected to emerge by 2030, nearly half of which is ecologically pristine (see figure 10). This pattern, in combination with factors encouraging rubber cultivation in higher altitudes and latitudes, underscores the urgency and importance of careful future land use planning. Local and regional decision-makers can use mid-to (more cautiously) long-term assessment such as this to develop policy guidelines and decision support mechanism that can take the occurrence of potential new land use and land management systems into account. Either to prepare a certain region for potential innovations regarding the demands to local infrastructure, or to put necessary guidelines and rules into place to 'soften the blow' these innovations might have on traditional systems or biodiversity and nature conservation.