# Forecasting the Influence of Climate Change on Agroecosystem Services: Potential Impacts on Honey Yields in a Small-Island Developing State

Authors: Diana L. Delgado; María Eglée Pérez; Alberto Galindo-Cardona; Tugrul Giray; Carla Restrepo
Journal: Psyche
 (2012)

---

## Abstract

Global change poses numerous challenges to developing nations and small-island developing states (SIDSs). Among these are the effects of climate change on honeybees’ provisioning services including honey production. Here we ask two questions. First, what is the relationship between honey yield and climate in a tropical environment? Second, how does yield vary spatially under current climate and future scenarios of climate change? Focusing on the island of Puerto Rico, we developed an ensemble of bioclimatic models that were used in a geographical information system to identify suitable areas for honey production under current and future scenarios of climate change. A comparison between contemporary (1998–2005) and historical (1910–1974) honey yield data revealed a reduction in average yield, including variability, over time, with current yields averaging 5.3 L/colony. Three bioclimatic variables were retained by at least three models: temperature seasonality and mean temperature of the wettest quarter were negatively correlated with honey yields whereas precipitation of the wettest month was positively correlated. The four models varied in terms of their predictions but showed that both honey yields and areas suitable for honey production will decrease under scenarios of climate change. These results illustrate the possible impacts of climate change on honey and ultimately honeybees.


---

## Body

## 1. Introduction

Honeybees (Apis mellifera) and their resource base have been managed to enhance supporting and provisioning services to human kind since ancient times [1–3]. In recent years much emphasis has been placed on the decline of pollination services provided by honeybees both in natural and managed ecosystems in response to multiple drivers of change [4]. This is because pollination directly impacts the functioning of ecosystems and ultimately local and regional economies. Surprisingly, little attention has been paid to potential changes in the delivery of provisioning services such as honey and beeswax production [5, 6]. Understanding these changes is important because beekeeping is promoted as a tool for rural development and conservation in developing nations in the tropics or regions therein (e.g., [7, 8]; http://www.beesfordevelopment.org/).The importance of honeybees in the delivery of provisioning services is reflected by the widespread introduction of beekeeping practices by Europeans to their colonies in the early 1600s [2]. In the Caribbean, beekeeping was aimed at the production of honey and beeswax sometimes to supply local, and at other times, regional, and international markets [9–11]. Yet the production of honey and beeswax seems to have varied greatly across the islands and within them as illustrated by the island of Puerto Rico where a “boom and bust” cycle occurred in tandem with a decline in honey yields (Figure 1(a); supplementary material available online at doi: 10.1155/2012/951215). Such cycles reflect the often complex dynamics of markets as influenced not only by the social, political, and technological vagaries of a region, but also the interactions between honeybees, their resource base, and environmental factors including meteorological events.(a) Historical (1910–1974) and (b) contemporary (1998–2005) honey yields for the island of Puerto Rico illustrating the “boom and bust” cycle of honey production. Based on historical accounts, it is possible to distinguish four periods in this cycle and the readers are referred to the supplementary material for details. Average (black diamond) and maximums (asterisks) honey yields, including standard deviations (boxes). In years for which honey yield data is sparse (≤8 municipalities) the standard deviation boxes are in white, otherwise in gray. The closed circles joined by a line depict the number of farms reporting honeybee colonies. Historical data were compiled from historical agricultural census data (supplementary material).
(a)
(b)Among the environmental factors that may impact the delivery of provisioning services by honeybees is climate change as the observed variation in honeybee abundance and honey yields along climatic gradients suggests [1, 12, 13]. At low latitudes, honeybees remain active throughout the year whereas at high latitudes they pass through a period of complete inactivity [1, 5, 14, 15]. Likewise within the tropics, the activity of honeybees decreases with increasing elevation [15]. Climate directly influences honeybee behavior given the strong dependency of bee foraging activity and flight on temperature, solar radiation, and wind at a variety of time scales [16, 17]. Climate can indirectly influence honeybees through its effects on their resource base, including flowering plants, pathogens, and predators [6, 18, 19]. Temperature and to a lesser degree precipitation seem to exert a primary control on honeybee activity, yet the extent to which climate change will impact honey yields is poorly understood.This lack of understanding of the effects of climate change on honey yields, and more broadly speaking the delivery of provisioning services, is prevalent atregional to local scales particularly in developing regions and small-island developing states (SIDSs) [20, 21]. First, there are large uncertainties regarding the effects of climate change scenarios at increasingly smaller scales due to the coarseness of climate change models and scarcity of climate data [22]. Second, ecosystem services are delivered at local scales but are influenced by processes operating at multiple scales [23]. Lastly, there is an uneven capacity among regions and nations to develop research and technology that could help cope adaptively with global change [24]. Developing regions and SIDS are a point in case. These regions not only face the greatest uncertainties [25–27] but also the greatest vulnerabilities [28] regarding crop production and trade under scenarios of global change. The general consensus is that the productivity of agroecosystems, in particular those delivering the major world food crops, will decrease towards the tropics and subtropics [29]. Therefore, it remains to be seen what happens with the vast majority of crops that are produced, consumed, and traded within regional and local markets, and at the same time are delivered by agroecosystems contributing additional ecosystem services. One such agroecosystem is honeybees and their resource base. Here we ask how honey yields will be impacted by climate change in the island of Puerto Rico, one SIDS within the Caribbean.We develop four models describing the relationship between honey yields and climate taking advantage of available contemporary (1998–2005) honey yield data for the island of Puerto Rico. Then we analyze these models spatially to identify the areas suitable for honey production under current conditions and future scenarios of climate change. Our approach is based on the development of an ensemble of bioclimatic forecasting models in which the combination of multiple forecasts increases the robustness of the predictions [30]. We reasoned that this approach was necessary given that several regionalized climate change models for the Caribbean show conflicting scenarios for the island of Puerto Rico [31–34]. Therefore, any attempt to evaluate the effects of climate change on the delivery of honey, and more broadly speaking ecosystem services in this region, needs to account for these diverse scenarios.


## 2. Methods

### 2.1. Study Site

Puerto Rico with its 8,740 km2 is the smallest island of the Greater Antilles in the Caribbean (Figure 2). The island has a diverse set of bioclimatic conditions due, in part, to its complex topography and wide elevation range (0–1,338 m; [35–39]). This diversity in combination with available data on contemporary beekeeping activity makes the island ideal for examining the relationship between honey yields and climate variability.Figure 2
Map of Puerto Rico showing its position within the Caribbean basin, the centroids of the municipalities with contemporary (1998–2005) honey yield data (circles) and the location of the weather stations (flags) with existing precipitation and maximum and minimum temperatures for the same period.The surplus of honey produced per colony or beehive, that is, honey yields, integrates management and environmental factors that directly and indirectly influence honeybee activity [40]. Modeling honey yields as a function of climate involved two steps: the compilation and creation of honey yield and climate spatial databases; the modeling of honey yields as a function of bioclimatic variables. The latter was generated in BIOCLIM, a predictive system developed for the purpose of modeling the distribution of animals and plants, including agriculturally important crops [41, 42].


### 2.2. Spatial Databases

We used unpublished beekeeping census data for the period 1998–2005 to calculate honey yields at the municipality level (Figure2). These data collected by the Department of Agriculture of Puerto Rico (DAPR), with some exceptions, mirror the historical agricultural census data (Figure 1(a)). The DAPR collects statistical data through surveys conducted every two years among beekeepers on the location of the beehives or bee colonies, the number of bee colonies, the total volume of honey produced, the amount of honey sold, and the income generated from the sale during the previous and the current year (A. M. Cruz pers. comm.). These unpublished census data exclude information on beeswax production. To maintain data confidentiality, DAPR does not disclose the exact location of bee colonies; instead they report the corresponding municipality. This may reduce the spatial resolution of the data and limit the possibility to verify unusual data points with the beekeepers. Nevertheless, this dataset provides valuable temporal and spatial information about honey yields on a regional scale (Figure 2).For each beekeeper, honey yields were calculated dividing the total volume of honey produced by the number of colonies that they reported. Then for each municipality we averaged honey yields across beekeepers to obtain yearly (1998–2005 period; the years 2001 and 2003 were excluded because of incomplete or missing data) and overall (6 years) honey yields. The mean number of beekeepers per municipality ranged between 2.4 and 2.8 and the average number of municipalities with honey yield data was30±7 out of a total of 78.We compiled monthly total precipitation and monthly average maximum and minimum temperatures from those weather stations whose records matched the time period covered by the DAPR’s honey yield datasets. A search of the National Oceanic and Atmospheric Administration’s (NOAA) cooperative weather stations (http://www.ncdc.noaa.gov/oa/ncdc.html), the Luquillo Long-Term Ecological Research Site (http://luq.lternet.edu/), and the Atmos Carib Research Center at the University of Puerto Rico-Mayaguez Campus (http://atmos.uprm.edu/) databases yielded 126 weather stations of which 21 met closely our criteria (Figure 2). Five of these 21 stations had continuous monthly data for the 6 years covered by this study (72 months). The remaining stations had ≤19 months with missing data and 2 among these had a full year of missing data (2003 and 2005 for Mayagüez city and Pico del Este stations). We completed the monthly missing data as follows: averaging monthly values across two adjacent years to complete months without data when daily data did not exist (Case 1), averaging daily values for a given month when incomplete daily data existed (Case 2), and predicting missing monthly data using linear regression models that related climate data from nearby stations (Case 3) in the case of Mayagüez city and Pico del Este stations.The honey yield and climatic data were added as attributes of point features in a GIS (State Plane Puerto Rico and Virgin Island FIPS 5200, NAD 83) to create interpolated surfaces with a 450 m resolution that were needed as input data in our spatially explicit modeling approach. In ArcGIS 9.2, we used the inverse distance-weighted (IDW) method with its default values (power = 2 and maximum number of neighbors = 15) to interpolate the yearly and overall honey yields added to the centroids of each municipality, the yearly and overall monthly total precipitation, and averaged monthly maximum and minimum temperatures added to each weather station (Figure2).We chose the IDW method over others (splice, kriging, natural neighbors) because it produced the interpolated surfaces that most resembled the actual climatic distribution in the island. IDW is a local, deterministic interpolation method that estimates unknown point values based on known neighboring sample points, whose influence decreases with distance according to a negative power function [43]. The error of the interpolated surfaces can vary with P and neighborhood characteristics and also with the variable under consideration [44–46]. We used a cross-validation procedure on some selected variables and varied P (1, 2, 3) and the maximum number of neighbors (10, 15) to examine their impact on the root mean square errors (RMSEs). We found that RMSE varied minimally, and therefore kept the interpolated surfaces created with the IDW default values.These interpolated climatic maps were used together with a digital elevation model (DEM; seamless.usgs.gov) in DIVA’s BIOCLIM module to generate yearly and overall bioclimatic variables (Table1; DIVA version 5.4.0.1; http://www.diva-gis.org/).Table 1
Bioclimatic variables calculated by BIOCLIM software and used in the stepwise multiple regressions to produce the four honey yield predictive models.
BIOCLIM variable
Description
Type
Bio 1
Annual mean temperature (°C)
Annual trend
Bio 21,2
Mean monthly range (°C)
Seasonality
Bio 31,2
Isothermality (Bio 2/Bio 7)
Seasonality
Bio 41,2,3,4
Temperature seasonality (CV×100)
Seasonality
Bio 51,2,3,4
Maximum temperature of warmest month (°C)
Extreme conditions
Bio 61,2,3,4
Minimum temperature of coldest month (°C)
Extreme conditions
Bio 7
Temperature annual range (Bio5-Bio6) (°C)
Seasonality
Bio 81,2,3,4
Mean temperature of wettest quarter (°C)
Extreme conditions
Bio 91,2,4
Mean temperature of driest quarter (°C)
Extreme conditions
Bio 10
Mean temperature of warmest quarter (°C)
Extreme conditions
Bio 11
Mean temperature of coldest quarter (°C)
Extreme conditions
Bio 12
Annual precipitation (mm)
Annual trend
Bio 131,2,3,4
Precipitation of wettest month (mm)
Extreme conditions
Bio 141,2,4
Precipitation of driest month (mm)
Extreme conditions
Bio 151,2,3,4
Precipitation seasonality (CV×100)
Seasonality
Bio 164
Precipitation of wettest quarter (mm)
Extreme conditions
Bio 171,2,3,4
Precipitation of driest quarter (mm)
Extreme conditions
Bio 181,2,3,4
Precipitation of warmest quarter (mm)
Extreme conditions
Bio 19
Precipitation of coldest quarter (mm)
Extreme conditions
A quarter represents a 3-month period.  The superscripts indicate the variables that were preselected by PCA and VIF analyses based on multicollinearity to enter the models (see text).


### 2.3. Modeling Contemporary Honey Yields

We developed an ensemble of four models that altogether provide a robust representation of current and future predicted honey yields [30]. The models reflected different ways of aggregating the data (yearly versus overall averages) and handling regions (all versus subset of municipalities) for which the bioclimatic data may not fully represent existing conditions. Model 1 included yearly data for all municipalities for which honeybee data was available. In model 2, we included yearly data for all but four municipalities located in eastern Puerto Rico, namely Ceiba, Naguabo, Luquillo, and Rio Grande. Although these municipalities encompass coastal areas that are characteristically dry, all operating weather stations in the area are clustered in the Luquillo Mountains that are characteristically wet. In Model 3, we averaged yearly values for all the municipalities for which we had honey yield data. Finally, in Model 4 we used the average yearly values of all variables as in Model 3 but excluded the same municipalities as in Model 2.To characterize each municipality based on honey yields and bioclimatic conditions, we averaged the pixel values of the corresponding raster maps and added these averages to their centroids (Figure2). For each modeled dataset we run exploratory data analysis (EDA) to help detect outliers or anomalies in the data. Subsequently, we ran a principal component analysis (PCA) on the correlation matrices both as a way to handle variables that were in different units and to preselect variables based on degree of multicollinearity [47]. To validate this approach, we calculated the variance inflation factor (VIF), a procedure that quantifies the extent of multicollinearity among multiple variables that will be included in regression models [47].The procedure outlined above helped us eliminate six variables. The remaining 13 variables were included in stepwise multiple regressions to explore the relationship between the log-transformed honey yield and bioclimatic variables (Table1). Both the Akaike’s information criterion (AIC) and the adjusted R-squares are used by Spotfire S+ (TIBCO) to select the model with the best fit. The AIC takes into account statistical goodness of fit while penalizing for increasing the number of variables; low AIC values indicate the model with the best fit [47]. Subsequently, we validated the models generated through the stepwise multiple regression analysis by using the automated model selection function dredgein the MuMIn package version 1.7.7 of R statistical software. This function examines all the possible models given the provided variables and ranks them according to their AIC [48].


### 2.4. Modeling Honey Yields under Scenarios of Climate Change

Various modeling efforts to predict climate change trends in the Caribbean agree in that the temperature will increase 2°C by 2099 but disagree regarding precipitation trends [31–34, 49]. For Puerto Rico, in particular, under the IS92 business as usual scenario (BaU), it was predicted that precipitation would not change during the dry season (DS; December–April) but that it would decrease 20 mm during the early rainy season (ERS; May–July) and increase 15 mm during the late rainy season (LRS; August–November) [32, 50]. In contrast, under the SRES A2 scenario, an ensemble of models predicted that the precipitation was going to decrease 28, 66, and 50 mm during the DS, ERS, and LRS seasons, respectively [31]. We used these two climate change scenarios to model honey yields and modified accordingly the original climate data used as input for BIOCLIM. The new maps depicting these conditions were used as input data to predict honey yields.We used ArcGIS 9.2 software to create the honey yields and climate spatial databases and to model honey yields across the island under current and future climate scenarios.


## 2.1. Study Site

Puerto Rico with its 8,740 km2 is the smallest island of the Greater Antilles in the Caribbean (Figure 2). The island has a diverse set of bioclimatic conditions due, in part, to its complex topography and wide elevation range (0–1,338 m; [35–39]). This diversity in combination with available data on contemporary beekeeping activity makes the island ideal for examining the relationship between honey yields and climate variability.Figure 2
Map of Puerto Rico showing its position within the Caribbean basin, the centroids of the municipalities with contemporary (1998–2005) honey yield data (circles) and the location of the weather stations (flags) with existing precipitation and maximum and minimum temperatures for the same period.The surplus of honey produced per colony or beehive, that is, honey yields, integrates management and environmental factors that directly and indirectly influence honeybee activity [40]. Modeling honey yields as a function of climate involved two steps: the compilation and creation of honey yield and climate spatial databases; the modeling of honey yields as a function of bioclimatic variables. The latter was generated in BIOCLIM, a predictive system developed for the purpose of modeling the distribution of animals and plants, including agriculturally important crops [41, 42].


## 2.2. Spatial Databases

We used unpublished beekeeping census data for the period 1998–2005 to calculate honey yields at the municipality level (Figure2). These data collected by the Department of Agriculture of Puerto Rico (DAPR), with some exceptions, mirror the historical agricultural census data (Figure 1(a)). The DAPR collects statistical data through surveys conducted every two years among beekeepers on the location of the beehives or bee colonies, the number of bee colonies, the total volume of honey produced, the amount of honey sold, and the income generated from the sale during the previous and the current year (A. M. Cruz pers. comm.). These unpublished census data exclude information on beeswax production. To maintain data confidentiality, DAPR does not disclose the exact location of bee colonies; instead they report the corresponding municipality. This may reduce the spatial resolution of the data and limit the possibility to verify unusual data points with the beekeepers. Nevertheless, this dataset provides valuable temporal and spatial information about honey yields on a regional scale (Figure 2).For each beekeeper, honey yields were calculated dividing the total volume of honey produced by the number of colonies that they reported. Then for each municipality we averaged honey yields across beekeepers to obtain yearly (1998–2005 period; the years 2001 and 2003 were excluded because of incomplete or missing data) and overall (6 years) honey yields. The mean number of beekeepers per municipality ranged between 2.4 and 2.8 and the average number of municipalities with honey yield data was30±7 out of a total of 78.We compiled monthly total precipitation and monthly average maximum and minimum temperatures from those weather stations whose records matched the time period covered by the DAPR’s honey yield datasets. A search of the National Oceanic and Atmospheric Administration’s (NOAA) cooperative weather stations (http://www.ncdc.noaa.gov/oa/ncdc.html), the Luquillo Long-Term Ecological Research Site (http://luq.lternet.edu/), and the Atmos Carib Research Center at the University of Puerto Rico-Mayaguez Campus (http://atmos.uprm.edu/) databases yielded 126 weather stations of which 21 met closely our criteria (Figure 2). Five of these 21 stations had continuous monthly data for the 6 years covered by this study (72 months). The remaining stations had ≤19 months with missing data and 2 among these had a full year of missing data (2003 and 2005 for Mayagüez city and Pico del Este stations). We completed the monthly missing data as follows: averaging monthly values across two adjacent years to complete months without data when daily data did not exist (Case 1), averaging daily values for a given month when incomplete daily data existed (Case 2), and predicting missing monthly data using linear regression models that related climate data from nearby stations (Case 3) in the case of Mayagüez city and Pico del Este stations.The honey yield and climatic data were added as attributes of point features in a GIS (State Plane Puerto Rico and Virgin Island FIPS 5200, NAD 83) to create interpolated surfaces with a 450 m resolution that were needed as input data in our spatially explicit modeling approach. In ArcGIS 9.2, we used the inverse distance-weighted (IDW) method with its default values (power = 2 and maximum number of neighbors = 15) to interpolate the yearly and overall honey yields added to the centroids of each municipality, the yearly and overall monthly total precipitation, and averaged monthly maximum and minimum temperatures added to each weather station (Figure2).We chose the IDW method over others (splice, kriging, natural neighbors) because it produced the interpolated surfaces that most resembled the actual climatic distribution in the island. IDW is a local, deterministic interpolation method that estimates unknown point values based on known neighboring sample points, whose influence decreases with distance according to a negative power function [43]. The error of the interpolated surfaces can vary with P and neighborhood characteristics and also with the variable under consideration [44–46]. We used a cross-validation procedure on some selected variables and varied P (1, 2, 3) and the maximum number of neighbors (10, 15) to examine their impact on the root mean square errors (RMSEs). We found that RMSE varied minimally, and therefore kept the interpolated surfaces created with the IDW default values.These interpolated climatic maps were used together with a digital elevation model (DEM; seamless.usgs.gov) in DIVA’s BIOCLIM module to generate yearly and overall bioclimatic variables (Table1; DIVA version 5.4.0.1; http://www.diva-gis.org/).Table 1
Bioclimatic variables calculated by BIOCLIM software and used in the stepwise multiple regressions to produce the four honey yield predictive models.
BIOCLIM variable
Description
Type
Bio 1
Annual mean temperature (°C)
Annual trend
Bio 21,2
Mean monthly range (°C)
Seasonality
Bio 31,2
Isothermality (Bio 2/Bio 7)
Seasonality
Bio 41,2,3,4
Temperature seasonality (CV×100)
Seasonality
Bio 51,2,3,4
Maximum temperature of warmest month (°C)
Extreme conditions
Bio 61,2,3,4
Minimum temperature of coldest month (°C)
Extreme conditions
Bio 7
Temperature annual range (Bio5-Bio6) (°C)
Seasonality
Bio 81,2,3,4
Mean temperature of wettest quarter (°C)
Extreme conditions
Bio 91,2,4
Mean temperature of driest quarter (°C)
Extreme conditions
Bio 10
Mean temperature of warmest quarter (°C)
Extreme conditions
Bio 11
Mean temperature of coldest quarter (°C)
Extreme conditions
Bio 12
Annual precipitation (mm)
Annual trend
Bio 131,2,3,4
Precipitation of wettest month (mm)
Extreme conditions
Bio 141,2,4
Precipitation of driest month (mm)
Extreme conditions
Bio 151,2,3,4
Precipitation seasonality (CV×100)
Seasonality
Bio 164
Precipitation of wettest quarter (mm)
Extreme conditions
Bio 171,2,3,4
Precipitation of driest quarter (mm)
Extreme conditions
Bio 181,2,3,4
Precipitation of warmest quarter (mm)
Extreme conditions
Bio 19
Precipitation of coldest quarter (mm)
Extreme conditions
A quarter represents a 3-month period.  The superscripts indicate the variables that were preselected by PCA and VIF analyses based on multicollinearity to enter the models (see text).


## 2.3. Modeling Contemporary Honey Yields

We developed an ensemble of four models that altogether provide a robust representation of current and future predicted honey yields [30]. The models reflected different ways of aggregating the data (yearly versus overall averages) and handling regions (all versus subset of municipalities) for which the bioclimatic data may not fully represent existing conditions. Model 1 included yearly data for all municipalities for which honeybee data was available. In model 2, we included yearly data for all but four municipalities located in eastern Puerto Rico, namely Ceiba, Naguabo, Luquillo, and Rio Grande. Although these municipalities encompass coastal areas that are characteristically dry, all operating weather stations in the area are clustered in the Luquillo Mountains that are characteristically wet. In Model 3, we averaged yearly values for all the municipalities for which we had honey yield data. Finally, in Model 4 we used the average yearly values of all variables as in Model 3 but excluded the same municipalities as in Model 2.To characterize each municipality based on honey yields and bioclimatic conditions, we averaged the pixel values of the corresponding raster maps and added these averages to their centroids (Figure2). For each modeled dataset we run exploratory data analysis (EDA) to help detect outliers or anomalies in the data. Subsequently, we ran a principal component analysis (PCA) on the correlation matrices both as a way to handle variables that were in different units and to preselect variables based on degree of multicollinearity [47]. To validate this approach, we calculated the variance inflation factor (VIF), a procedure that quantifies the extent of multicollinearity among multiple variables that will be included in regression models [47].The procedure outlined above helped us eliminate six variables. The remaining 13 variables were included in stepwise multiple regressions to explore the relationship between the log-transformed honey yield and bioclimatic variables (Table1). Both the Akaike’s information criterion (AIC) and the adjusted R-squares are used by Spotfire S+ (TIBCO) to select the model with the best fit. The AIC takes into account statistical goodness of fit while penalizing for increasing the number of variables; low AIC values indicate the model with the best fit [47]. Subsequently, we validated the models generated through the stepwise multiple regression analysis by using the automated model selection function dredgein the MuMIn package version 1.7.7 of R statistical software. This function examines all the possible models given the provided variables and ranks them according to their AIC [48].


## 2.4. Modeling Honey Yields under Scenarios of Climate Change

Various modeling efforts to predict climate change trends in the Caribbean agree in that the temperature will increase 2°C by 2099 but disagree regarding precipitation trends [31–34, 49]. For Puerto Rico, in particular, under the IS92 business as usual scenario (BaU), it was predicted that precipitation would not change during the dry season (DS; December–April) but that it would decrease 20 mm during the early rainy season (ERS; May–July) and increase 15 mm during the late rainy season (LRS; August–November) [32, 50]. In contrast, under the SRES A2 scenario, an ensemble of models predicted that the precipitation was going to decrease 28, 66, and 50 mm during the DS, ERS, and LRS seasons, respectively [31]. We used these two climate change scenarios to model honey yields and modified accordingly the original climate data used as input for BIOCLIM. The new maps depicting these conditions were used as input data to predict honey yields.We used ArcGIS 9.2 software to create the honey yields and climate spatial databases and to model honey yields across the island under current and future climate scenarios.


## 3. Results

### 3.1. Honey Yields

Average honey yields for the period 1998–2005 were estimated at5.3±4.4 L/colony (mean ± SD) whereas for the historical data this figure was 11.3±9.6 L/colony (t-test, df = 465, P≤0.05; Figure 1). Contemporary honey yields seem to be less variable (smaller standard deviations (shown above) and have lower maximums (20.5 in 2002 versus 78.7 L/colony in 1910)) than the historical values (Figure 1). Honey yields for the period 1998–2005 as well as the historical dataset exhibit a large intra- and inter-annual variability that may reflect to a large extent differences among the municipalities and time periods in terms of their socioeconomic and ecological potential to sustain honeybee activity. We use the average honey yield (5.3 L/colony) for the 1998–2005 period as a baseline figure to compare the behavior of the models that predict honey yields under current and future scenarios of climate change.


### 3.2. Modeling Current Honey Yields

The four models differed in terms of the total variance that they explained and the bioclimatic variables associated with honey yields (Table2). Models 1 and 2 explained the least amount of variance in the data (R2~0.18) whereas the opposite was true for Models 3 and 4 (R2~0.50). The models also differed regarding the type and number of variables that were retained with Model 1 and Model 3 retaining 2 and 8, respectively. Four models retained temperature seasonality (Bio 4) and mean temperature of the wettest quarter (Bio 8), three retained precipitation of the wettest month (Bio 13), and two retained the minimum temperature of the coldest month (Bio 6) (Table 2). Honey yields were negatively correlated with temperature seasonality (Bio 4) and mean temperature of the wettest quarter (Bio 8), and positively correlated with minimum temperature of the coldest month (Bio 6) and the precipitation of the wettest month (Bio 13).Table 2
Results of stepwise multiple regressions with predictors retained in each of the four predictive models of honey production.
Model
Model parameter
1(all years and municipalities)
2(all years but not all municipalities)
3(yearly average and all municipalities)
4(yearly average but not all municipalities)
Intercept
3.6328
4.6138
5.1130
6.3048
Isothermality (Bio 3)
—
−0.0099
—
—
Temperature seasonality (Bio 4)
−0.0021
−0.0043
−0.0037
−0.0083
Maximum temperature of warmest month (Bio 5)
—
—
0.1130
—
Minimum temperature of coldest month (Bio 6)
—
—
0.0715
0.0594
Mean temperature of wettest quarter (Bio 8)
−0.0647
−0.0688
−0.2382
−0.1863
Precipitation of wettest month (Bio 13)
—
0.0004
0.0039
0.0012
Precipitation seasonality (Bio 15)
—
—
−0.0285
—
Precipitation of driest quarter (Bio 17)
—
—
−0.0033
—
Precipitation of warmest quarter (Bio 18)
—
—
−0.0009
—
P
<0.001
<0.001
<0.001
<0.001
R
2
0.1903
0.1847
0.5276
0.4984Predicted current honey yields varied among the models (Figures3 and 4, and Table 3). In Models 1 and 3, these ranged between 1.0 and 14.0 L/colony whereas in Model 2 and 4 between 2.0 and 67.0 L/colony, more than doubling the maximum yields predicted by the former (Table 3). The location and extent of the areas suitable for honey yields ≥5.3 L/colony varied among the four models. Models 3 and 4 identified areas suitable for honey production not shown by the other two, and as a result the degree of overlap among the four models was 35%; eliminating Model 4, the most dissimilar model, gave an overlay of 48%. Finally, the predicted areal extent of areas suitable for honey production varied among the models and ranged between 1,000 and 2,200 km2 or equivalently between 11 and 25% of the total area of the island (Model 1 < Model 3 < Model 2 < Model 4; Table 3).Table 3
Predicted honey yields and areal extents of suitable honey production areas (honey yields ≥ 5.3 L/colony).  Under each model the future climate change scenarios are based on (A) Neelin et al. [31], which follows the A2 scenario, and (B) Angeles et al. [32], the IS92 business as usual scenario.
Model
Predicted honey yields (L/colony)
Predicted areal extents of suitable honey production areas (km2)
Current
Future
A
B
Current
Future
Max–min (Mean ± SD)
Max–min (Mean ± SD)
Max–min (Mean ± SD)
A
B
1
3.5–14.0(4.6±0.9)
2.6–10.6(3.5±0.7)
2.6–10.7(3.6±0.7)
1,012
221
250
2
3.2–32.3(5.1±2.1)
2.3–23.7(3.7±1.6)
2.4–24.2(3.9±1.5)
1,538
617
704
3
0.6–12.0(4.4±1.4)
0.2–10.5(3.4±1.4)
0.6–11.6(4.1±1.3)
1,370
917
1,149
4
1.7–67.0(6.2±6.2)
0.9–39.1(3.5±3.8)
1.0–39.4(3.9±3.7)
2,210
782
909Figure 3
Shade relief map of Puerto Rico showing the predicted current areas suitable for honey yields ≥5.3 L/colony.Figure 4
Maps showing the calculated estimates of current honey yield (models 1, 2, 3, and 4) and the honey yields predicted under future climate change scenarios. The future scenarios include a 2°C increase in temperature with a 28, 66, and 50 mm decrease in precipitation during the DS, ERS, and LRS, respectively (Future A) and a 2°C increase in temperature with no change in DS, 20 mm decrease in precipitation during ERS, and 15 mm increase in precipitation during LRS (Future B).


### 3.3. Modeling Honey Yields under Climate-Change Scenarios

Under the two scenarios of climate change, the predicted honey yields ranged between 10.5 and 39.4 L/colony, which are lower than the predicted current values (Table3). Yet under the SRES’s A2 scenario the minimum and maximum honey yields were slightly lower than under the ISP2’s BaU scenario, a result that becomes obvious when examining the areas suitable for honey production (Figure 4; Table 3). All models predicted a reduction in the areal extent of areas suitable for honey production. In the SRES’s A2 scenario this reduction ranges between 33% and 78% and under the ISP2’s BaU scenario between 16% and 75% (Models 3 < 2 < 4 < 1).


## 3.1. Honey Yields

Average honey yields for the period 1998–2005 were estimated at5.3±4.4 L/colony (mean ± SD) whereas for the historical data this figure was 11.3±9.6 L/colony (t-test, df = 465, P≤0.05; Figure 1). Contemporary honey yields seem to be less variable (smaller standard deviations (shown above) and have lower maximums (20.5 in 2002 versus 78.7 L/colony in 1910)) than the historical values (Figure 1). Honey yields for the period 1998–2005 as well as the historical dataset exhibit a large intra- and inter-annual variability that may reflect to a large extent differences among the municipalities and time periods in terms of their socioeconomic and ecological potential to sustain honeybee activity. We use the average honey yield (5.3 L/colony) for the 1998–2005 period as a baseline figure to compare the behavior of the models that predict honey yields under current and future scenarios of climate change.


## 3.2. Modeling Current Honey Yields

The four models differed in terms of the total variance that they explained and the bioclimatic variables associated with honey yields (Table2). Models 1 and 2 explained the least amount of variance in the data (R2~0.18) whereas the opposite was true for Models 3 and 4 (R2~0.50). The models also differed regarding the type and number of variables that were retained with Model 1 and Model 3 retaining 2 and 8, respectively. Four models retained temperature seasonality (Bio 4) and mean temperature of the wettest quarter (Bio 8), three retained precipitation of the wettest month (Bio 13), and two retained the minimum temperature of the coldest month (Bio 6) (Table 2). Honey yields were negatively correlated with temperature seasonality (Bio 4) and mean temperature of the wettest quarter (Bio 8), and positively correlated with minimum temperature of the coldest month (Bio 6) and the precipitation of the wettest month (Bio 13).Table 2
Results of stepwise multiple regressions with predictors retained in each of the four predictive models of honey production.
Model
Model parameter
1(all years and municipalities)
2(all years but not all municipalities)
3(yearly average and all municipalities)
4(yearly average but not all municipalities)
Intercept
3.6328
4.6138
5.1130
6.3048
Isothermality (Bio 3)
—
−0.0099
—
—
Temperature seasonality (Bio 4)
−0.0021
−0.0043
−0.0037
−0.0083
Maximum temperature of warmest month (Bio 5)
—
—
0.1130
—
Minimum temperature of coldest month (Bio 6)
—
—
0.0715
0.0594
Mean temperature of wettest quarter (Bio 8)
−0.0647
−0.0688
−0.2382
−0.1863
Precipitation of wettest month (Bio 13)
—
0.0004
0.0039
0.0012
Precipitation seasonality (Bio 15)
—
—
−0.0285
—
Precipitation of driest quarter (Bio 17)
—
—
−0.0033
—
Precipitation of warmest quarter (Bio 18)
—
—
−0.0009
—
P
<0.001
<0.001
<0.001
<0.001
R
2
0.1903
0.1847
0.5276
0.4984Predicted current honey yields varied among the models (Figures3 and 4, and Table 3). In Models 1 and 3, these ranged between 1.0 and 14.0 L/colony whereas in Model 2 and 4 between 2.0 and 67.0 L/colony, more than doubling the maximum yields predicted by the former (Table 3). The location and extent of the areas suitable for honey yields ≥5.3 L/colony varied among the four models. Models 3 and 4 identified areas suitable for honey production not shown by the other two, and as a result the degree of overlap among the four models was 35%; eliminating Model 4, the most dissimilar model, gave an overlay of 48%. Finally, the predicted areal extent of areas suitable for honey production varied among the models and ranged between 1,000 and 2,200 km2 or equivalently between 11 and 25% of the total area of the island (Model 1 < Model 3 < Model 2 < Model 4; Table 3).Table 3
Predicted honey yields and areal extents of suitable honey production areas (honey yields ≥ 5.3 L/colony).  Under each model the future climate change scenarios are based on (A) Neelin et al. [31], which follows the A2 scenario, and (B) Angeles et al. [32], the IS92 business as usual scenario.
Model
Predicted honey yields (L/colony)
Predicted areal extents of suitable honey production areas (km2)
Current
Future
A
B
Current
Future
Max–min (Mean ± SD)
Max–min (Mean ± SD)
Max–min (Mean ± SD)
A
B
1
3.5–14.0(4.6±0.9)
2.6–10.6(3.5±0.7)
2.6–10.7(3.6±0.7)
1,012
221
250
2
3.2–32.3(5.1±2.1)
2.3–23.7(3.7±1.6)
2.4–24.2(3.9±1.5)
1,538
617
704
3
0.6–12.0(4.4±1.4)
0.2–10.5(3.4±1.4)
0.6–11.6(4.1±1.3)
1,370
917
1,149
4
1.7–67.0(6.2±6.2)
0.9–39.1(3.5±3.8)
1.0–39.4(3.9±3.7)
2,210
782
909Figure 3
Shade relief map of Puerto Rico showing the predicted current areas suitable for honey yields ≥5.3 L/colony.Figure 4
Maps showing the calculated estimates of current honey yield (models 1, 2, 3, and 4) and the honey yields predicted under future climate change scenarios. The future scenarios include a 2°C increase in temperature with a 28, 66, and 50 mm decrease in precipitation during the DS, ERS, and LRS, respectively (Future A) and a 2°C increase in temperature with no change in DS, 20 mm decrease in precipitation during ERS, and 15 mm increase in precipitation during LRS (Future B).


## 3.3. Modeling Honey Yields under Climate-Change Scenarios

Under the two scenarios of climate change, the predicted honey yields ranged between 10.5 and 39.4 L/colony, which are lower than the predicted current values (Table3). Yet under the SRES’s A2 scenario the minimum and maximum honey yields were slightly lower than under the ISP2’s BaU scenario, a result that becomes obvious when examining the areas suitable for honey production (Figure 4; Table 3). All models predicted a reduction in the areal extent of areas suitable for honey production. In the SRES’s A2 scenario this reduction ranges between 33% and 78% and under the ISP2’s BaU scenario between 16% and 75% (Models 3 < 2 < 4 < 1).


## 4. Discussion

Average annual honey yields in the island of Puerto Rico were estimated at 5.3 L/colony. The four models developed to predict honey yields under current and future scenarios of climate change varied in terms of the predicted honey yields and the extent and location of areas suitable for honey production (area suitable for honey production; yields ≥5.3 L/colony). The predicted current (minimum and maximum ranges were 0.6–3.5 and 12–67.0 L/colony, resp.) and future (0.2–2.6 and 10.5–39.4 L/colony, resp.) honey yields indicate that climate change has the potential to reduce yields almost by half. In addition, the predicted areal extents of the current (minimum and maximum range between 1,000 and 2,200 km2) and future (221–1,149 km2) areas suitable for honey production show a substantial decrease further supporting the likely impact of climate change on beekeeping. Overall these results indicate that climate change has the potential to affect the delivery of provisioning and supporting services by honeybees.


### 4.1. Honey Yields

Current average honey yields are almost half the historical ones (5.3 versus 11.3 L/colony; Figure1) and lower than the world (10.7 and 13.1 L/colony for the years 1984 and 1998, resp.) and Caribbean (17.3 and 14.6 L/colony, resp.) estimates ([1]; Food Agricultural Organization (FAO) 1998 database, http://faostat.fao.org/). The long-term trends in honey yields and number of farms with beehives reconstructed from a variety of historical sources resemble a “boom and bust” cycle in which complex interactions between humans, honeybees, and their resource base determine fluctuations in the delivery of ecosystem services (Figure 1; supplementary material). These long-term trends also resemble a population growth model that overshoots its carrying capacity and crashes. In Puerto Rico, a reduction in maintenance research characteristic of modern agriculture [51], in combination with multiple socioeconomic [52, 53] and environmental factors including food availability, diseases, invasions, and natural meteorological events, such as hurricanes, explain these long-term trends and contemporary low honey yields in the island (supplementary material).Understanding these long-term trends is important because it raises questions about cycles of production and the magnitude of climate change impacts depending on the stage along these cycles. Furthermore, it raises questions about the characteristics of these cycles in developing regions and SIDS (Figure1).


### 4.2. Modeling Current Honey Yields

We identified a subset of bioclimatic variables that explains part of the variability in current honey yields, as well as its spatial variability in a tropical region. Three temperature- and one precipitation-derived variables were common to ≥2 models, highlighting the influence of temperature on honey yields. Specifically, honey yields decreased with greater temperature seasonality (Bio 4) and mean temperature of the wettest quarter (Bio 8) in all four models and increased with precipitation of the wettest month (Bio 13) and the minimum temperature of the coldest month (Bio 6) in three and two models, respectively. Although there is considerable variability in these and other bioclimatic variables across the island [35, 37], the general trend is for the wettest quarter of the year that runs from July to September to be a period with the highest maximum temperatures of the year. Likewise, the coldest month is January, a month with some of the lowest precipitations, and the wettest month October, a month in which the high summer temperatures begin to ease. Our results suggest that a combination of extreme temperatures and low precipitation limits the activity of honeybees in this tropical setting.The map of predicted honey yields identifies regions in eastern Puerto Rico and in the Central Mountains as suitable for average or above average honey yields, which further clarifies the interpretation of our results. These mountainous regions have milder temperatures and are more humid than the lowlands, in particular the dry, low-laying areas of southern and western Puerto Rico [35–39]. Most likely then, the observed variability in predicted honey yields reflects ways in which different life history components of honeybees vary along the complex bioclimatic gradients observed in the island. In Costa Rica, for example, honeybees remain active throughout the year, yet worker brood increased and honey production decreased along an elevation gradient (900–2800 m) located in a mesic environment [15]. This effect, however, was more pronounced during the rainy than the dry season. In Germany, the weight of bee colonies varied across a vast region in Hannover largely in response to the observed variability in climatic conditions and land use [14].


### 4.3. Modeling Honey Yields under Scenarios of Climate-Change

Honey yields as well as the area suitable for honey production decreased under both climate change scenarios. These results are not surprising given the strong dependency of honeybees on temperature even in tropical regions [5, 14, 15]. More interesting were the effects of precipitation on predicted honey yields. In each of the four models, the minimum and maximum yields did not differ markedly between the two climate change scenarios, but the same was not true for the predicted area suitable for honey production (honey yields ≥5.3 L/colony). In a warmer and drier Puerto Rico, the area will decrease whereas in a warmer and wetter Island the change will be lessened. In Model 3, the model in which several precipitation variables were retained, the area deemed suitable for production was the largest observed. In the model 3b, where there is an overall precipitation increase, this effect was more pronounced. One possibility is that a precipitation increase converts some of the dry areas of the island into areas suitable for honey production (Table 3, Figure 4).Results of our work mirror those obtained by others focusing on major food crops delivered by agroecosystems in developing regions and SIDS [25–27]. Specifically it has been shown a decrease in yields towards the tropics [54]. Given the long-term and often complex dynamics of agricultural production, we may ask if climate change in interaction with socioeconomic and environmental factors will contribute to further yield declines [55]. Based on Puerto Rico’s beekeeping history we can speculate that the extent to which climate change can impact this and similar agroecosystems including the services that we derive from them may depend upon the stage at which they are found along the “boom and bust” cycle.


### 4.4. Future Work

We only included bioclimatic variables in our modeling efforts, and the fact that the models explained upto 53% of the variability in honey yields suggests that other variables should be considered. Topography, wind speed, colony management, and land use are logical factors to explore. The first two may influence honeybee behavior as already outlined in our Introduction whereas the third might play a minor role in Puerto Rico’s current setting for reasons already mentioned. Land use on the other hand, may integrate characteristics of the food resource base utilized by honeybees, including its quality, quantity, and spatial distribution [6, 14]. It would be worth exploring more mechanistic models that could directly assess the impact of raising temperature and atmospheric CO2 levels on the physiology of honeybees and their food sources. In particular, raising atmospheric CO2 levels may impact this agroecosystem in different ways. First, it may increase plant photosynthetic activity and water use efficiency [56], thus increasing the bees’ food source base. Second, in some flowering species increases in CO2 may cause a reduction in nectar production [57], the principal source of food for bees. Lastly, exposure of honeybees to high levels of CO2 for prolong periods of time can alter insect physiology and behavior, and in some cases it becomes lethal [58, 59].


### 4.5. Application and Relevance

Beekeeping is promoted as a tool for rural development and conservation not only in the Caribbean but other tropical nations or regions therein ([7, 8]; http://www.beesfordevelopment.org/). Likewise, honeybees are valued for the pollination services provided both to agricultural [4, 60–62] and natural ecosystems, including endangered species [63, 64]. Therefore, any initiative that may increase people’s dependency on honeybees in developing regions and SIDS should take into account the likely effect of climate change on beekeeping. Already there are regions where drought is mentioned as the main problem for beekeeping because it leads to shortages of bees’ food followed by colony absconding [65]. Equally, shortages of food may increase negative interactions between native bees and honeybees due to niche overlap. Our work is an example of an approach that can provide a better understanding of the bioclimatic factors that limit honey production, and by doing so it may help farmers to cope with new environmental conditions.


## 4.1. Honey Yields

Current average honey yields are almost half the historical ones (5.3 versus 11.3 L/colony; Figure1) and lower than the world (10.7 and 13.1 L/colony for the years 1984 and 1998, resp.) and Caribbean (17.3 and 14.6 L/colony, resp.) estimates ([1]; Food Agricultural Organization (FAO) 1998 database, http://faostat.fao.org/). The long-term trends in honey yields and number of farms with beehives reconstructed from a variety of historical sources resemble a “boom and bust” cycle in which complex interactions between humans, honeybees, and their resource base determine fluctuations in the delivery of ecosystem services (Figure 1; supplementary material). These long-term trends also resemble a population growth model that overshoots its carrying capacity and crashes. In Puerto Rico, a reduction in maintenance research characteristic of modern agriculture [51], in combination with multiple socioeconomic [52, 53] and environmental factors including food availability, diseases, invasions, and natural meteorological events, such as hurricanes, explain these long-term trends and contemporary low honey yields in the island (supplementary material).Understanding these long-term trends is important because it raises questions about cycles of production and the magnitude of climate change impacts depending on the stage along these cycles. Furthermore, it raises questions about the characteristics of these cycles in developing regions and SIDS (Figure1).


## 4.2. Modeling Current Honey Yields

We identified a subset of bioclimatic variables that explains part of the variability in current honey yields, as well as its spatial variability in a tropical region. Three temperature- and one precipitation-derived variables were common to ≥2 models, highlighting the influence of temperature on honey yields. Specifically, honey yields decreased with greater temperature seasonality (Bio 4) and mean temperature of the wettest quarter (Bio 8) in all four models and increased with precipitation of the wettest month (Bio 13) and the minimum temperature of the coldest month (Bio 6) in three and two models, respectively. Although there is considerable variability in these and other bioclimatic variables across the island [35, 37], the general trend is for the wettest quarter of the year that runs from July to September to be a period with the highest maximum temperatures of the year. Likewise, the coldest month is January, a month with some of the lowest precipitations, and the wettest month October, a month in which the high summer temperatures begin to ease. Our results suggest that a combination of extreme temperatures and low precipitation limits the activity of honeybees in this tropical setting.The map of predicted honey yields identifies regions in eastern Puerto Rico and in the Central Mountains as suitable for average or above average honey yields, which further clarifies the interpretation of our results. These mountainous regions have milder temperatures and are more humid than the lowlands, in particular the dry, low-laying areas of southern and western Puerto Rico [35–39]. Most likely then, the observed variability in predicted honey yields reflects ways in which different life history components of honeybees vary along the complex bioclimatic gradients observed in the island. In Costa Rica, for example, honeybees remain active throughout the year, yet worker brood increased and honey production decreased along an elevation gradient (900–2800 m) located in a mesic environment [15]. This effect, however, was more pronounced during the rainy than the dry season. In Germany, the weight of bee colonies varied across a vast region in Hannover largely in response to the observed variability in climatic conditions and land use [14].


## 4.3. Modeling Honey Yields under Scenarios of Climate-Change

Honey yields as well as the area suitable for honey production decreased under both climate change scenarios. These results are not surprising given the strong dependency of honeybees on temperature even in tropical regions [5, 14, 15]. More interesting were the effects of precipitation on predicted honey yields. In each of the four models, the minimum and maximum yields did not differ markedly between the two climate change scenarios, but the same was not true for the predicted area suitable for honey production (honey yields ≥5.3 L/colony). In a warmer and drier Puerto Rico, the area will decrease whereas in a warmer and wetter Island the change will be lessened. In Model 3, the model in which several precipitation variables were retained, the area deemed suitable for production was the largest observed. In the model 3b, where there is an overall precipitation increase, this effect was more pronounced. One possibility is that a precipitation increase converts some of the dry areas of the island into areas suitable for honey production (Table 3, Figure 4).Results of our work mirror those obtained by others focusing on major food crops delivered by agroecosystems in developing regions and SIDS [25–27]. Specifically it has been shown a decrease in yields towards the tropics [54]. Given the long-term and often complex dynamics of agricultural production, we may ask if climate change in interaction with socioeconomic and environmental factors will contribute to further yield declines [55]. Based on Puerto Rico’s beekeeping history we can speculate that the extent to which climate change can impact this and similar agroecosystems including the services that we derive from them may depend upon the stage at which they are found along the “boom and bust” cycle.


## 4.4. Future Work

We only included bioclimatic variables in our modeling efforts, and the fact that the models explained upto 53% of the variability in honey yields suggests that other variables should be considered. Topography, wind speed, colony management, and land use are logical factors to explore. The first two may influence honeybee behavior as already outlined in our Introduction whereas the third might play a minor role in Puerto Rico’s current setting for reasons already mentioned. Land use on the other hand, may integrate characteristics of the food resource base utilized by honeybees, including its quality, quantity, and spatial distribution [6, 14]. It would be worth exploring more mechanistic models that could directly assess the impact of raising temperature and atmospheric CO2 levels on the physiology of honeybees and their food sources. In particular, raising atmospheric CO2 levels may impact this agroecosystem in different ways. First, it may increase plant photosynthetic activity and water use efficiency [56], thus increasing the bees’ food source base. Second, in some flowering species increases in CO2 may cause a reduction in nectar production [57], the principal source of food for bees. Lastly, exposure of honeybees to high levels of CO2 for prolong periods of time can alter insect physiology and behavior, and in some cases it becomes lethal [58, 59].


## 4.5. Application and Relevance

Beekeeping is promoted as a tool for rural development and conservation not only in the Caribbean but other tropical nations or regions therein ([7, 8]; http://www.beesfordevelopment.org/). Likewise, honeybees are valued for the pollination services provided both to agricultural [4, 60–62] and natural ecosystems, including endangered species [63, 64]. Therefore, any initiative that may increase people’s dependency on honeybees in developing regions and SIDS should take into account the likely effect of climate change on beekeeping. Already there are regions where drought is mentioned as the main problem for beekeeping because it leads to shortages of bees’ food followed by colony absconding [65]. Equally, shortages of food may increase negative interactions between native bees and honeybees due to niche overlap. Our work is an example of an approach that can provide a better understanding of the bioclimatic factors that limit honey production, and by doing so it may help farmers to cope with new environmental conditions.


## 5. Conclusions

Contemporary and historical beekeeping records in Puerto Rico revealed a likely reduction in average honey yields, including variability, over time. Current honey yields were used as baseline to compare the behavior of an ensemble of four spatially explicit models that were developed to predict honey yields under current and future scenarios of climate change. The four models varied in their predictions, yet they all showed that honey yields, as well as the area suitable for honey production, will decrease under scenarios of climate change. These results illustrate the possible impacts of climate change on honeybees and ultimately the essential services that they provide to us.



---

*Source: 951215-2012-10-30.xml*