temperature, precipitation), which are combined to determine areas and environments that adhere to the species niche, which are manifest as positive areas on the map. The BIOCLIM variables [49] (bioclimatic variables representing annual trends, seasonality and extremes of the environment; accumulated from weather station data from the 1961–1990), accessed via the WorldClim website (http://www.worldclim.org/bioclim. Accessed 2012 Jan 25), were used as grid data at 30 arc seconds resolution (approximately 1 km61 km). To reduce the processing overhead and data storage, the data was clipped to Africa only. We used the default settings in MaxEnt, which have been shown to give robust and reliable results [44]. The initial models of predicted distribution were tested by omitting 20% of the samples (39 sample points for testing and 158 to produce the model), to gauge the robustness of the models. Area under the curve (AUC), jackknife, analysis of variant contribution, and response curves were also calculated by MaxEnt for further analysis. Preliminary iterations of the models allowed us to identify doubtful outlying localities (i.e. where predictions were very low), and these were corrected where necessary (e.g. geo-referencing errors). A low predictive value was identified for the populations on Mt. Marsabit (Kenya), and so two independent data sets were used, one with and one without Mt. Marsabit. Spatial sorting bias can give a false impression of model fitness by inflating AUC values [50]. We initially accounted for spatial sorting bias by removing all duplicate localities, (i.e. those with identical georeferences: 713 data records were reduced to 349 unique localities, representing a reduction of 48.9% of the sampling), then these 349 localities were reduced to 197 samples (only one per 1 km grid cell) for the MaxEnt modelling. Then, to more fully test for spatial sorting bias and to check the robustness of the models, further models were produced with spatial sorting bias removed, and for only the core region. Spatial sorting bias was removed by only including points that were separated by at least 0.2 of a degree (c. 22 km), and then checked using ArcGIS 10 [41] average nearest neighbour, which gives a clustering qualification, from clustered, to random, to dispersed. The original 197 samples were highly cluster, but after re-sampling (22 samples points remaining) the points were spatially random. These 22 points were remodelled in MaxEnt, both for the whole of Africa and for a very restricted core region. It should be noted that removal of the bias in this way was only done to test the fitness (comparing AUC and results) of the models; the rest of the analysis was performed using the 349 unique localities (of which 197 were used as the samples for MaxEnt). Predictive maps represent potential maps of distribution, not the actual distribution of a species: they show where the ecological niche is potentially suitable for the species. The species may not actually occur in this predicted area, for any of the following reasons. (1) The predicted area will include niches that are no longer suitable for the species, because the vegetation has been drastically altered (e.g. deforestation, agriculture). (2) There is niche saturation, e.g. another species or many species already occupy that niche. (3) The models are incorrect due to inaccuracies and weak resolution in the recorded environmental data, which is often the case for microclimates. It should be made clear that the original climate data is also modelled, and based on data from a network of global weather stations. (4) The layers used for modelling are not necessarily the drivers for the occurrence of the species, and other environmental factors (which may drive occurrence) have not been included in the model (e.g. soil characteristics, at a suitable resolution). (5) The locality data is inaccurate. For indigenous Arabica, at least the last two of these points should not strongly influence our distribution model, Climate Change and Indigenous Arabica Coffee because the predictive model is very robust (see Results) and the data have been carefully verified by careful scrutiny using satellite maps and ground-truthing, respectively. It should also be said that Arabica coffee is closely tied to narrow environmental parameters, and like the vast majority of coffee species [22] it has a restricted and specific distribution. Data for cultivated Arabica confirms these observations, and shows that this species is sensitive to environmental variables, particularly temperature and precipita- tion [2,6,8]. Soil-water balance is a key factor for the growth of Arabica, although until now this data is only available for the species in cultivation [2,8]. In addition, observation and reported distribution show that coffee species are rarely adventive, are only scarcely found in secondary (regenerating) vegetation, and quickly become stressed in degraded habitats [22]. In general, the family Rubiaceae includes a very high number of narrowly endemic species, which in most species appears to be related to ecological sensitivity [51]. Future mapping/climate change modelling Future predictive maps for indigenous Arabica were generated using the same BIOCLIM layers as employed in the MaxEnt analysis, with predicted future climate data downscaled using the delta method to 30 arc seconds (c. 1 km61 km). The delta method interpolates the General Circulation Model generally used in climate modelling at scales of 100 to 200 km using a thin plate spline spatial interpolation method to achieve the 30 arc seconds resolution [52]. The data was provided by the CGIAR Research Program on Climate Change, Agriculture and Food Security (CCAFS) (http:// www.ccafs-climate.org/statistical downscaling_delta/) [52]. The Met Office Hadley Centre (Hadley Centre for Climate Prediction and Research) climate change model, Hadley Centre Coupled Model, version 3 (HadCM3) [53], a coupled atmosphere-ocean general circulation model, was used for the time intervals 2020, 2050 and 2080 (note these date represent a time windows of ten years either side of the time interval date, i.e. 2020 is an average of the years 2010–2029, 2050 for 2040–2059 and 2080 for 2070–2089), under three emission scenarios of the IPCC Special Report on Emissions Scenarios (SRES) [54]: scenario A1B (maximum energy requirements; emissions differentiated dependent on fuel sources; balance across sources), A2A (high energy requirements; emissions less than A1/Fl) and B2A (lower energy requirements; emissions greater than B1). We tried, and subsequently rejected, the consensus forecasting with unweighted average approach [55], because it did not provide any further quality to the analyses. There are several climate models available (e.g. HadCM3, CCCMA and CSIRO) but the HadCM3 model was selected for this study because it is presently the only one that spans the years 2020, 2050 and 2080, for all three scenarios, at the resolution of 1 km. It is also reported to provide good median results for Africa compared with other models [13]. Separate model iterations were run with and without Mt Marsabit (Kenya), nine for each iteration, giving a total of 18 models. Locality data (i.e. the 349 unique localities) were queried against all model combinations, in order to provide a value for each of the localities for each of the models. For the present day models, based on all localities, thresholds were calculated to cut off at 68%, 95% and 100% of the 349 (unique) localities, i.e. based on their relative values within the niche model. As we do not have a totally normal distribution for the prediction at each sample point, 68% is not quite equal to 1 standard deviation. These thresholds were used to classify each Arabica locality though time, within each scenario, into 68% (optimal), 95% (intermediate [suboptimal]; includes the 68% threshold), 100% (marginal [extreme environments]; includes 68% and 95% thresholds). These thresholds were applied to the models PLOS ONE | www.plosone.org 3 November 2012 | Volume 7 | Issue 11 | e47981