Climate change impact on coffee and the pollinator bee suitable area interaction in Raya Azebo, Ethiopia

Modeling suitable habitat of plants and pollinators helps to find plant and pollinator interactions. Flowereng plants are dependent on insect pollinators, pollintors are also dependent on flowering plants. However, there is mismatch of plant and pollinator suitable area due to deficiency of plentiful research. Therefore, this study investigates the suitable area of coffee, bee, and their interaction. This was conducted using bee and coffee presence, altitude, and climate variables. This study was conducted using spatial regression and ecological niche model. In addition, future common suitable areas were detected in whole Ethiopia using climate analogue in ensemble general circulation models (GCMs) in order to adapt the impact of climate change. Thus, the result showed that there is 666 km suitable area for coffee production. Precipitation seasonality and precipitation of warmest quarter are the most contributors for coffee distribution. Similarly, a suitable area for bee is 835 km. Precipitation of driest quarter and precipitation of warmest quarter are the best contributors for bee distribution. Their common suitable area in relation to the enviromental variables is 598 km. Their interaction is explained by 52.93%. The model shows that there might be future common suitable area in Ethiopia. Therefore, the study remarks that coffee and bee should be cultivated in their common suitable area in order to ensure their pollination, ecological and economic roles under changing climate. Subjects: Environment & Agriculture; Bioscience; Environmental Studies & Management


Introduction
Agricultural and ecological change affects smallholder farmers (Vermeulen et al., 2012). Plant and animal distributions are shifting to find suitable environments as climate changes (Thomas ABOUT THE AUTHOR Haftu Abrha is a lecturer at Institute of Climate Society (ICS) in Mekelle University. His area of interest is in forestry, envoronment, species distribution modeling and climate related studies.

PUBLIC INTEREST STATEMENT
This article draws on the experience of potential suitable area of coffee and bee in relation to environmental variables and their interactions in order to indicate plant-pollinator interaction. Plant-pollinator interactions can be interrupted through spatial mismatches of plant and pollinator distribution. Therefore, this article highlights the need for development of plant and pollinator interaction studies to improve ecosystem services of insect pollinators and to conserve biodiversity under changing climate. et al., 2004). Climate change affects the geographic distribution and effectiveness of pollinators (Giannini et al., 2012).
Pollination is an essential input in crop production compared to any other input of crop production such as fertilizer, weeding, labor, or pesticides (Klein et al., 2006). Most crops produce in the presence of pollinators (Ricketts, 2004). Insect pollinators are important for transporting genes of plant species (Kearns, Inouye, & Waser, 1998). Among the insects, bees afford indispensable pollination services to many flowering plants (Wratten, Gillespie, Decourtye, Mader, & Desneux, 2012). Beekeeping activity can be done in agroforestry and marginal areas. However, early identification of beekeeping suitable areas is an important activity to increase productivity and long-term survival of the colony (Malczewski, 2006). Recently, much attention focused on the management of honeybees because their decline is a serious threat to biodiversity (Zoccali et al., 2017). Managed honeybees can provide important ecosystem services by increasing pollination facility to agriculture and wild plants (Cox-Foster et al., 2007).
Ethiopia has good potential for honey production because of its favorable conditions for beekeeping. Ethiopia produces around 23.58% and 2.13% of the total share of honey production on African and on a global, respectively (Abebe & Puskur, 2011). Beekeeping enhances productivity and biodiversity (Ajao & Oladimeji, 2013). However, drought, pests and disease, shortage of bee forage, and lack of skill were the major constraints in Ethiopia (Abebe & Puskur, 2011). In addition, climate change has an impact on beekeeping (Langowska et al., 2017). Climate change declines honeybees and causes reducing of pollinator activity (Hegland, Nielsen, Lázaro, Bjerknes, & Totland, 2009). Therefore, bees are under threat due to climate change and environmental degradation (Potts et al., 2016).
Coffee is the world's most prized tropical export crop (Craparo, Van Asten, Läderach, Jassogne, & Grab, 2015). It is grown in more than 60 tropical countries on over 11 million hectares by an estimated 25 million farmers (Waller, Bigger, & Hillocks, 2007). Ethiopia is the center of origin of coffee and is the largest country producing Arabica coffee (Wondimkun, Jebessa, Molloro, & Haile, 2016). Coffee is contributing the lion's share in the national economy because it is a source of foreign exchange earnings (Kufa, Ayano, Yilma, Kumela, & Tefera, 2011). The major coffeeproducing regions in the country are Oromiya and Southern Nations Nationalities and Peoples Region, whereas other regions like Tigray produce coffee in small amounts (Sualeh & Mekonnen, 2013). In Tigray, Raya Azebo, coffee is grown with khat in farm level.
Agroforestry systems, commonly used to produce coffee, provide an important ecosystem service such as increasing biodiversity and productivity (Jha et al., 2014). Recent studies predict severe climate change impacts on coffee production and distribution (Craparo et al., 2015). Coffee is affected by climate change in two ways: directly, through extreme climate events on coffee production and indirectly, through changes in pollination efficiency (Cox-Foster et al., 2007).
Plant-pollinator interactions can be interrupted through spatial mismatches of plant and pollinator distribution. Climate change affects the phenology, abundance, and distribution of plants and pollinators. There is still limited knowledge of how climate change affects plant-pollinator mutualisms and how changed availability of mutualistic partners influences the persistence of interacting species (Hegland et al., 2009). Therefore, suitability analyses can be extremely helpful to plan land uses, according to specific requirements (Malczewski, 2006). This is used to identify the potential area of coffee-bee suitable areas using species distribution models (SDMs) and spatial regression models. SDMs are used to model current and future suitable areas of pollinators (Giannini et al., 2017) and plants (Abrha, Birhane, Hagos, & Manaye, 2018). Hence, the study aims to model potential suitable area of coffee and honeybee in relation to environmental variables, investigate the spatial interaction of coffee and honeybee in relation to environmental variables, and to identify future common suitable areas of coffee and bee under changing climate using climate analog tool.

Study area description
Raya Azebo is located between 12°17ʹ and 12°15ʹ north and 39°22ʹ to 39°25ʹ east ( Figure 1). The mean annual temperature of the district is 21.5°C. Its mean annual rainfall is 779 mm. Its rainfall is a bimodal type pattern from February to April period and from July to September. There is high potential evapotranspiration (Bewket, Radeny, & Mungai, 2015). The agro-climatic condition of the area is characterized as midland (85%) and lowland (15%) (Lemlem, 2013).

Species distribution modeling
SDMs are currently the most broadly used scientific methods to study potential climate change impacts on biodiversity (Franklin et al., 2013). The models combine species geographical coordinates with a set of predictor variables. The models produces a species distribution map using the physical characteristics of the study area (Miola, Freitas, Barbosa, & Fernandes, 2011). Species distribution models (SDMs) are displaying good results (Ceccarelli et al., 2015).
In recent times, several types of models have been used to simulate the potential distribution of species. These models included BIOCLIM, Random Forest, Boosted Regression Tree, BIOMAPPER, DIVA, and DOMAIN, CLIMEX, GAM, GLM, GARP, and MaxEnt (Kriticos & Randall, 2001;Shabani, Kumar, & Ahmadi, 2016). Of these, MaxEnt was selected because of the advantages of using presence-only data and performing well with incomplete data, small sample sizes, and gaps . Maximum entropy algorithm (Maxent) is one of the most accurate and globally used ecological niche model (Phillips, Anderson, & Schapire, 2006). Modeling potential species distributions using Maxent requires the use of environmental variables to recognize any changes in their distribution (Evangelista, Kumar, Stohlgren, & Young, 2011). Models with more than 30 location points are capable of yielding accurate and precise results (Baldwin, 2009). Therefore, 75 bee and 38 coffee latitude/longitude occurrence points were collected with the aid of a handheld GPS.

Bioclimatic variables
World climate (WordClim) is a set of global climate layers with gridded climate datasets (from 1950 to 2000). It consists of 19 global gridded data fields created from average monthly temperature and precipitation values. It is with a spatial resolution from 30 s (~1 km 2 ) to 10 min (~340 km 2 ). It is available in world climate website (http://www.worldclim.org/tiles.php?Zone=27). To produce a detailed and high-resolution image, climatic variables with 30 arc second resolution were used. The database includes 19 bioclimatic variables (Table 1).
A large amount of studies are conducted using bioclimatic variables to model areas suitable for species across spatial and temporal scales (Carnaval, Hickerson, Haddad, Rodrigues, & Moritz, 2009;Carstens & Richards, 2007;Nogués-Bravo, 2009) especially in environmental, agricultural, and biological sciences (Parra, Graham, & Freile, 2004). Current climate data and altitude downloaded from the word climate database were clipped down in GIS to the Ethiopia-Raya Azebo map extension to model the current distribution.

Model performance measurement
Occurrence points were divided randomly into training data and test data to test the predictive performance of SDMs. Of these points, 25% were used to validate the model, and the remaining 75% were used to calibrate the model. The accuracy of the model was measured by area under the curve (AUC) of the receiver operating characteristic (ROC) plot. AUC value is between 0.5 and 1.0 to determine the accuracy of the model (Mbatudde, Mwanjololo, Kakudidi, & Dalitz, 2012).

Environmental variable contribution
Environmental variable contribution on the species was assessed through percentage contribution table and jackknife test (Phillips, Anderson, & Schapire, 2006). The percent contribution of each variable to the final model involves a series of trials until the model gain increases (Baldwin, 2009). As well, jackknife measures the importance of each variable by examining the influence on the prediction, in isolation (Girma, de Bie, Skidmore, Venus, & Bongers, 2016). Finally, the relationship between the potential distribution and environmental factors is known.

Suitability threshold
The Maxent model output is a continuous raster of probability values ranging between 0 and 1 which represents habitat suitability. The Maxent algorithm produces results ranging from 0 to 1, demonstrating relative suitability of a given grid (Graham & Hijmans, 2006, Phillips et al., 2006. Accordingly, a decision was made to categorize the distribution as either suitable or unsuitable (Evangelista et al., 2011). An area with p < 0.1772 (unsuitable area), 0.1772 ≤ p < 0.3543 (less suitable area), 0.3543 ≤ p < 0.5315 (suitable area), 0.5315 ≤ p < 0.7087 (optimum area), and p ≥ 0.7087 (excellent area), Reddy et al. (2014). These threshold values were classified in Arc and Diva GIS. 0.1772 was used to describe unsuitable and less to excellent suitable area.

Regression analysis
Regression analyses attempt to demonstrate the degree to which one or more variables potentially promote positive or negative change in another variable. Regression analysis is used to investigate and model the relationship between variables (Uyanık & Güler,). Diva GIS was used to calculate simple and multiple spatial regression. Simple regression was used to model coffee and bee interaction. Multiple regression was used to develop a model between common area of coffee and bee with other 20 environmental variables. Multiple regression gives insight into the relative importance of each environmental variable on the species presence (Peeters & Gardeniers, 1998). This was used to identify the contribution and significant of environmental variables on a common area of bee-coffee interaction. The level of significance of the variables was tested at a less strict 95% confidence level. In addition, F-statistic and critical values were used to test the level of significance between variables.

Climate analogues
Climate analogue technique detects areas where the climate today is similar to the future projected climate of another area (Luedeling & Neufeldt, 2012). This compares locations based on similarity in climate and other input variables (soils, social, and economic conditions). Climate analogue can predict through backward (where can I find the future climate of my site today?) and forward (where can I find the present climate of my site in the future?). This is used to adapt the impact of climate change on a species.
The climate analogue tool announced through Climate Change And Food Security (CCAFS) was used to identify analogous climates across space and time. The CCAFS model was employed to spatially locate analogue sites, i.e., a common suitable area of coffee-bee interaction of Raya Azebo based on the reliability (Diffenbaugh & Giorgi, 2012). The model was built using location points of the study area, 30 second resolution (1 km 2 ), 24 ensemble GCMs in a forward direction, and one greenhouse gas emissions scenarios (a1b) for the period of 2020 to 2049 (Hallegatte, Hourcade, & Ambrosi, 2007). A1b scenario is intermediate emission scenario which is similar to representative concentration pathways (RCP) 6 from the IPCC's Special Report on Emission Scenarios (Moss et al., 2010, Van Vuuren et al., 2011. Mean temperature and precipitation, projected from the CCAFS analog tiles, were the variables used in this study. Each of these variables is assigned an equal weight of 0.5 to their respective contributions. A 50% probability threshold was used to cut off low analogue site predictions, and then the raster having high similarity was selected for future distribution.

Results and discussions
4.1. Current distribution of coffee and bee

Model performance
The average test AUC value for the replicate runs of bee suitable area was 0.996 with standard deviation of 0.000. In addition, the average test AUC value for the replicate runs of coffee suitable area was 0.998 with standard deviation of 0.001 (Figure 2). A model with the highest AUC value is considered as the best performer (Sutherst, 2014). Pearce and Ferrier (2000) stated that models with an AUC >0.7 are applicable. These values indicated that the constructed model had 'Excellent' predictive accuracy.

Analysis of determinant variables
Precipitation of warmest quarter, precipitation of driest quarter, temperature seasonality, annual precipitation, and precipitation of driest month are top five environmental predictors intended for bee distribution. The environmental variable with the highest gain when used in isolation is precipitation of warmest quarter, which therefore appears to have the most useful information by itself. The environmental variable that decreases the gain the most when it is omitted is precipitation of warmest quarter, which therefore appears to have the most information that is not present in the other variables.
Precipitation seasonality, precipitation of warmest quarter, precipitation of coldest quarter, temperature annual range, and precipitation of driest month are the top five environmental predictors for coffee distribution. The environmental variable with the highest gain when used in isolation is precipitation of driest quarter, which therefore appears to have the most useful information by itself. The environmental variable that decreases the gain the most when it is omitted is precipitation of driest quarter, which therefore appears to have the most information that is not present in the other variables (Table 2).
A jackknife test points out the contribution of environmental variables similar to the percent contribution. The jackknife test is illustrated by a bar graph of dark blue, red, and lighter blue bars which indicates the level of contribution of environmental variables on the species distribution ( Figure 3).

Coffee and bee distribution
The potential distribution of the species contains areas of high and low predicted habitat suitability. The suitable sitesof bee and coffee are predicted in the western and northern part of Raya Azebo district. Besides, a suitable area of bee is in the eastern and northern parts of the study area. The current distribution of coffee is pointed out in the green color which covers about 666 km 2 (4a). In addition, the current distribution of bee is located in the red-yellowish color which covers about 835 km 2 (Figure 4(b)).
The study finds out the combined and individual suitable area of coffee and bee for agricultural planning processes. The lack of information about spatiotemporal change is common uncertainty in agricultural planning processes (Littleboy, Smith, & Bryant, 1996). Detecting suitable habitats to decide land use kind can be used to manage these uncertainties.
Similar to this study, Ovalle-Rivera, Läderach, Bunn, Obersteiner, and Schroth (2015) used coffee presence data and climate information data to predict climatic suitability for coffee production. Estrada, Rasche, & Schneider (2017) finds that climate variables have a high influence on coffee suitability. Accordingly, climate variability reduces the extent of suitable area (Hannah et al., 2017). Similarly, it shifts to new suitable areas.  There is a decrease of abundance of pollinators significantly due to lack of identification of suitable area (Ricketts, 2004). Therefore, this study identifies suitable areas of bee because insect pollinators increase the production, such as the number of pods, seeds per pod and seed weights per plant (Atmowidi, Buchori, Manuwoto, Suryobroto, & Hidayat, 2007), and fruit and seed sets of plants (Rianti, Suryobroto, & Atmowidi, 2010).

Coffee and bee common suitable area
The common suitable area for coffee and bee is depicted in green color in the map and covers about 598 km 2 . It is found in the northwestern part of the study area ( Figure 5(a)). The relationship between coffee and bee occurrence is significantly correlated with a coefficient of determination (R 2 ) of 0.5293. The p value is <0.00001, and result is significant at p < 0.05. This indicates that 52.93% of bee presence is explained by coffee presence. The remaining 47.07% is explained by other variables. Intercept and slope between coffee and bee are 0.0753 and 0.791, respectively ( Figure 5(b)). In addition, F-statistic value (F 1, 2078 = 2337.4) greater than the critical value (3.85) indicates that the relationship is significant. Similarly, a coefficient of determination is used to measure how well the regression predictions approximate the real data points (Asuero, Sayago, & Gonzalez, 2006). Hence, this indicates the model was good.
This study investigated that there is a significant relationship between bee and coffee. Similar to this study, bees have a strong interaction with flowering plants (Hanna, Foote, & Kremen, 2013). Insect pollinators are affected due to by the presence/absence of flowering plant and climatic conditions (Kajobe & Echazarreta, 2005). Arenas and Farina (2012) investigated the significant mutual relationships between insect pollinators and flowering plants. Mutualistic biotic interactions between flowering plants and their insect pollinators are a key component of biodiversity (Stein et al., 2017). This  interaction is fundamental for success of both pollinator and plant species (Landry, 2012). Insects, primarily bees, are responsible for most of the pollination services for wild and domestic plants (Potts et al., 2010). However, there are temporal and spatial mismatches between plants and their pollinators because of limited number of studies (Ashworth, Quesada, Casas, Aguilar, & Oyama, 2009).

Contribution of environmental variable on coffee and bee common area
The environmental variables having a significant contribution on coffee and bee common suitable area are mean temperature of wettest quarter, mean temperature of driest quarter, annual precipitation, precipitation of warmest quarter, precipitation of coldest quarter, and altitude (Table 3).

Conclusion and recommendations
In light of the above results and discussions, it is concluded that precipitation seasonality and precipitation of warmest quarter are the greatest contributors for coffee distribution. Precipitation of driest quarter and precipitation of warmest quarter are the best contributors for bee distribution. Mean temperature of wettest quarter, mean temperature of driest quarter, annual precipitation, precipitation of warmest quarter, precipitation of coldest quarter, and altitude have a significant contribution on coffee-bee common suitable area. Coffee and bee have significant interaction. The model shows there might be a future suitable area for both species.
Therefore, the study recommends that coffee and bee should be cultivated in their common suitable area in order to ensure their long-term survival, pollination, and ecological and economic roles. The study has several limitations in the datasets including soil, vegetation, and topography. Accordingly, additional datasets might improve the study.

Funding
The author received no direct funding for this research.

Competing Interests
The authors declare no competing interests.