Psyche

Volume 2012 (2012), Article ID 951215, 10 pages

doi:10.1155/2012/951215

Research Article

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

Diana L. Delgado,¹ María Eglée Pérez,² Alberto Galindo-Cardona,¹ Tugrul Giray,¹ and Carla Restrepo¹

¹Department of Biology, University of Puerto Rico-Rio Piedras Campus, P.O. Box 23360, San Juan, PR 00931-3360, USA

²Department of Mathematics, University of Puerto Rico-Rio Piedras Campus, P.O. Box 23355, San Juan, PR 00931-3355, USA

Received 16 May 2012; Accepted 19 August 2012

Academic Editor: Zachary Huang

Copyright © 2012 Diana L. Delgado et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

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.

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)

(b)

Figure 1: (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).

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 at regional 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 km² 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 (Figure 2). 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 was 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 (Figure 2).

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 and neighborhood characteristics and also with the variable under consideration [44–46]. We used a cross-validation procedure on some selected variables and varied (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 (Table 1; 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 ()
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 ()
Seasonality

Bio 16⁴
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 (Figure 2). 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 (Table 1). Both the Akaike’s information criterion (AIC) and the adjusted -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 dredge in the MuMIn package version 1.7.7 of 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 at  L/colony (mean ± SD) whereas for the historical data this figure was  L/colony (-test, df = 465, ; 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 (Table 2). Models 1 and 2 explained the least amount of variance in the data () whereas the opposite was true for Models 3 and 4 (). 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.4984

Predicted current honey yields varied among the models (Figures 3 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 km² 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 (km²)

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
2.6–10.6
2.6–10.7

221
250

2
3.2–32.3
2.3–23.7
2.4–24.2

617
704

3
0.6–12.0
0.2–10.5 
0.6–11.6

917

4
1.7–67.0
0.9–39.1
1.0–39.4

782
909

Figure 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 (Table 3). 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 km²) and future (221–1,149 km²) 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; Figure 1) 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 (Figure 1).

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 CO₂ levels on the physiology of honeybees and their food sources. In particular, raising atmospheric CO₂ 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 CO₂ may cause a reduction in nectar production [57], the principal source of food for bees. Lastly, exposure of honeybees to high levels of CO₂ 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.

Authors’ Contribution

D. L. Delgado Conceived and designed the study, performed the research, analyzed the data, and wrote the paper; M. E. Pérez Analyzed the data; A. Galindo-Cardona conceived or designed the study; T. Giray contributed new methods; C. Restrepo conceived and designed the study, performed the research, analyzed the data, and wrote the paper.

Acknowledgments

The authors give special thanks to Ana María Cruz and Edna Negrón Gorgas of the Department of Agriculture of Puerto Rico for providing the beekeeping census data, to D. Neelin and J. Gonzalez for providing climate scenarios data for Puerto Rico, Reinaldo Caraballo for sharing information on beekeeping in Puerto Rico, Nivea Santiago for locating historical articles on beekeeping, and Johanna Colón, Andrew Kerkhoff, Kimberley With, Jess Zimmerman, and several anonymous reviewers for providing useful comments on the paper. This project was partially funded through NSF-HRD-CREST (CR) and an NSF-PRLSAMP graduate fellowship (DD).

References

E. Crane, Bees and Beekeeping. Science, Practice, and World Resources, Comstock Publishing Associates, Ithaca, NY, USA, 1990.

E. Crane, The World History of Beekeeping and Honey Hunting, Routledge, New York, NY, USA, 1999.

C. Kremen, N. M. Williams, M. A. Aizen et al., “Pollination and other ecosystem services produced by mobile organisms: a conceptual framework for the effects of land-use change,” Ecology Letters, vol. 10, no. 4, pp. 299–314, 2007.

D. W. Roubik, “Tropical agriculture: the value of bees to the coffee harvest,” Nature, vol. 417, no. 6890, article 708, 2002.

W. Holmes, “The influence of weather on annual yields of honey,” Journal of Agricultural Science, vol. 139, no. 1, pp. 95–102, 2002.

X. Janssens, É. Bruneau, and P. Lebrun, “Prévision des potentialités de production de miel à l’échelle d’un rucher au moyen d’un système d’information géographique,” Apidologie, vol. 37, no. 3, pp. 351–365, 2006.

J. C. H. Saha, “Beekeeping for rural development, its potentiality and beekeeping against poverty-Bangladesh perspective,” in Proceedings of the 38th Congress Apimondia, Apimondia, Ljubljana, Slovenia, 2003.

J. B. Shrestha, “Combating desertification with sericulture and apiculture,” in Agriculture and Environment, pp. 6–13, Ministry of Agriculture and Cooperatives, Kathmandu, Nepal, 2006.

R. De La Sagra, Historia Económica-Política y Estadística de la Isla de Cuba o sea de sus Progresos en la Población, la Agricultura, el Comercio, y las Rentas, Imprenta de las Viudas de Arazoza y Soler, La Habana, Cuba, 1831.

C. M. Trelles, Biblioteca Científica Cubana, Imprenta de J. F. Olivier, Matanzas, Cuba, 1918.

A. I. Root Company, Cuba as a Bee Country, A. I. Root Company, Medina, Ohio, USA, 1903.

D. van Engelsdorp, J. Hayes, R. M. Underwood, and J. Pettis, “A survey of honey bee colony losses in the U.S., Fall 2007 to Spring 2008,” PLoS ONE, vol. 3, no. 12, Article ID e4071, 2008.

Y. L. Conte and M. Navajas, “Climate change: impact on honey bee populations and diseases,” OIE Revue Scientifique et Technique, vol. 27, no. 2, pp. 485–510, 2008.

J. Gerlach, “The dependence of honey-bee behaviour on weather—shown by the results of weighing bee-hives in the period of 1969 to 1978 in the region of the Beekeepers Association of Hannover,” International Journal of Biometeorology, vol. 29, no. 1, pp. 67–85, 1985.

M. Spivak, “The relative success of Africanized and European honey-bees over a range of life-zones in Costa Rica,” Journal of Applied Ecology, vol. 29, no. 1, pp. 150–162, 1992.

N. Vicens and J. Bosch, “Weather-dependent pollinator activity in an apple orchard, with special reference to Osmia cornuta and Apis mellifera (Hymenoptera: Megachilidae and Apidae),” Environmental Entomology, vol. 29, no. 3, pp. 413–420, 2000.

Z. Puškadija, E. Štefanić, A. Mijić et al., “Influence of weather conditions on honey bee visits (Apis mellifera carnica) during sunflower (Helianthus annuus L.) blooming period,” Poljoprivreda, vol. 13, no. 1, pp. 230–233, 2007.

N. L. Bradley, A. C. Leopold, J. Ross, and W. Huffaker, “Phenological changes reflect climate change in Wisconsin,” Proceedings of the National Academy of Sciences of the United States of America, vol. 96, no. 17, pp. 9701–9704, 1999.

A. M. Murilhas, “Varroa destructor infestation impact on Apis mellifera carnica capped worker brood production, bee population and honey storage in a Mediterranean climate,” Apidologie, vol. 33, no. 3, pp. 271–281, 2002.

S. R. Carpenter, R. DeFries, T. Dietz et al., “Millennium ecosystem assessment: research needs,” Science, vol. 314, no. 5797, pp. 257–258, 2006.

Millenium Ecosystem Assessment, Ecosystems and Human Well-Being: Synthesis, World Resources Institute, Washington, DC, USA, 2005.

UNFCCC, Climate Change: Impacts, Vulnerabilities, and Adaptation in Developing Countries, United Nations Framework Convention on Climate Change:, Bonn, Germany, 2007.

D. Schröter, W. Cramer, R. Leemans et al., “Ecology: ecosystem service supply and vulnerability to global change in Europe,” Science, vol. 310, no. 5752, pp. 1333–1337, 2005.

R. S. J. Tol, T. E. Downing, O. J. Kuik, and J. B. Smith, “Distributional aspects of climate change impacts,” Global Environmental Change, vol. 14, no. 3, pp. 259–272, 2004.

J. Keane, et al., Climate Change and Developing Country Agriculture: An Overview of Expected Impacts, Adaptation and Mitigation Challenges, and Funding Requirements, ICTSD-IPC Platform on Climate Change, Agriculture and Trade, Issue Brief No. 2, International Centre for Trade and Sustainable Development, and International Food & Agricultural Trade Policy Council, Washington, DC, USA, 2009.

P. G. Jones and P. K. Thornton, “The potential impacts of climate change on maize production in Africa and Latin America in 2055,” Global Environmental Change, vol. 13, no. 1, pp. 51–59, 2003.

P. K. Thornton, P. G. Jones, G. Alagarswamy, and J. Andresen, “Spatial variation of crop yield response to climate change in East Africa,” Global Environmental Change, vol. 19, no. 1, pp. 54–65, 2009.

FAO, Small Island Developing States. Agricultural Production and Trade, Preferenes, and Policy, Food and Agriculture Organization of the United Nations:, Rome, Italy, 2005.

M. L. Parry, C. Rosenzweig, A. Iglesias, M. Livermore, and G. Fischer, “Effects of climate change on global food production under SRES emissions and socio-economic scenarios,” Global Environmental Change, vol. 14, no. 1, pp. 53–67, 2004.

M. B. Araújo and M. New, “Ensemble forecasting of species distributions,” Trends in Ecology and Evolution, vol. 22, no. 1, pp. 42–47, 2007.

J. D. Neelin, M. Münnich, H. Su, J. E. Meyerson, and C. E. Holloway, “Tropical drying trends in global warming models and observations,” Proceedings of the National Academy of Sciences of the United States of America, vol. 103, no. 16, pp. 6110–6115, 2006.

M. E. Angeles, J. E. Gonzalez, D. J. Erickson, and J. L. Hernández, “Predictions of future climate change in the caribbean region using global general circulation models,” International Journal of Climatology, vol. 27, no. 5, pp. 555–569, 2007.

J. D. Campbell, et al., “Future climate of the Caribbean from a regional climate model,” International Journal of Climatology, vol. 31, no. 12, pp. 1866–1878, 2010.

ECLAC, Regional Climate Modeling in the Caribbean, Economic Commission for Latin America and the Caribbean-ECLAC, United Nations, 2010.

O. L. Fassig, “Average annual rainfall of Porto Rico, W. I,” Monthly Weather Review, vol. 37, pp. 982–986, 1909.

C. Daly, E. H. Helmer, and M. Quinones, “Mapping the climate of Puerto Rico, Vieques and Culebra,” International Journal of Climatology, vol. 23, no. 11, pp. 1359–1381, 2003.

B. A. Malmgren and A. Winter, “Climate zonation in Puerto Rico based on principal components analysis and an artificial neural network,” Journal of Climate, vol. 12, no. 4, pp. 977–985, 1999.

M. Jury, B. A. Malmgren, and A. Winter, “Subregional precipitation climate of the Caribbean and relationships with ENSO and NAO,” Journal of Geophysical Research D, vol. 112, no. 16, Article ID D16107, 2007.

J. J. Ewel and J. L. Whitmore, “The ecological life zones of Puerto Rico and the U.S. virgin islands,” Forest Service Research Paper ITF-18, 1973.

E. E. Southwick and D. Pimentel, “Energy efficiency of honey production by bees,” Bioscience, vol. 31, pp. 730–732, 1981.

A. Villordon, W. Njuguna, S. Gichuki et al., “Using GIS-Based tools and distribution modeling to determine sweetpotato germplasm exploration and documentation priorities in sub-Saharan Africa,” HortScience, vol. 41, no. 6, pp. 1377–1381, 2006.

L. J. Beaumont, L. Hughes, and M. Poulsen, “Predicting species distributions: use of climatic parameters in BIOCLIM and its impact on predictions of species' current and future distributions,” Ecological Modelling, vol. 186, no. 2, pp. 250–269, 2005.

A. D. Hartkamp, et al., Interpolation Techniques for Climate Variables, NRG-GIS Series 99-01, CIMMYT, Texcoco, Mexico, 1999.

A. Basistha, D. S. Arya, and N. K. Goel, “Spatial distribution of rainfall in Indian Himalayas—a case study of Uttarakhand Region,” Water Resources Management, vol. 22, no. 10, pp. 1325–1346, 2008.

D. Kurtzman and R. Kadmon, “Mapping of temperature variables in Israel: a comparison of different interpolation methods,” Climate Research, vol. 13, no. 1, pp. 33–43, 1999.

I. A. Nalder and R. W. Wein, “Spatial interpolation of climatic Normals: test of a new method in the Canadian boreal forest,” Agricultural and Forest Meteorology, vol. 92, no. 4, pp. 211–225, 1998.

G. P. Quinn and M. J. Keough, Experimental Design and Data Analysis for Biologists, Cambridge University Press, Cambridge, UK, 2004.

K. Barton, “Multi-model inference,” 2012

A. Chen, et al., “Climate trends and scenarios for climate change in the insular Caribbean,” Tech. Rep. 381, Caribbean Natural Resources Institute, Port Spain, Trinidad and Tobago, 2008.

M. E. Angeles, J. E. Gonzalez, D. J. Erickson, and J. L. Hernández, “An assessment of future Caribbean climate changes using the BAU scenario by coupling a global circulation model with a regional model,” in Proceedings of the 86th AMS Annual Meeting, Atlanta, Ga, USA, February 2006.

D. L. Plucknett and N. J. H. Smith, “Sustaining agricultural yields,” Bioscience, vol. 36, pp. 40–45, 1986.

M. A. Aizen and L. D. Harder, “Geographic variation in the growth of domesticated honey bee stocks: disease or economics?” Communitative and Integrative Biology, vol. 2, no. 6, pp. 464–466, 2009.

K. M. Roddy and L. Arita-Tsutsumi, “A history of honey bees in the Hawaiian Islands,” Journal of Hawaiian and Pacific Agriculture, vol. 8, pp. 59–70, 1997.

C. Rosenzweig and M. L. Parry, “Potential impact of climate change on world food supply,” Nature, vol. 367, no. 6459, pp. 133–138, 1994.

G. Fischer, M. Shah, F. N. Tubiello, and H. Van Velhuizen, “Socio-economic and climate change impacts on agriculture: an integrated assessment, 1990–2080,” Philosophical Transactions of the Royal Society B, vol. 360, no. 1463, pp. 2067–2083, 2005.

B. G. Drake, M. A. Gonzàlez-Meler, and S. P. Long, “More efficient plants: a consequence of rising atmospheric CO₂?” Annual Review of Plant Biology, vol. 48, pp. 609–639, 1997.

H. P. Rusterholz and A. Erhardt, “Effects of elevated CO₂ on flowering phenology and nectar production of nectar plants important for butterflies of calcareous grasslands,” Oecologia, vol. 113, no. 3, pp. 341–349, 1998.

G. Nicolas and D. Sillans, “Immediate and latent effects of carbon dioxide on insects,” Annual Review of Entomology, vol. 34, pp. 97–116, 1989.

X. Wang, L. Ji, G. Wang, and Y. Liu, “Potential effects of elevated carbon dioxide on leaf-feeding forest insects,” Frontiers of Biology in China, vol. 3, no. 1, pp. 68–77, 2008.

B. M. Freitas and R. J. Paxton, “The role of wind and insects in cashew (Anacardium occidentale) pollination in NE Brazil,” Journal of Agricultural Science, vol. 126, no. 3, pp. 319–326, 1996.

T. H. Ricketts, G. C. Daily, P. R. Ehrlich, and C. D. Michener, “Economic value of tropical forest to coffee production,” Proceedings of the National Academy of Sciences of the United States of America, vol. 101, no. 34, pp. 12579–12582, 2004.

D. W. Roubik, Pollination of Cultivated Plants in the Tropics, Food and Agriculture Organization of the United Nations, Rome, Italy, 1995.

M. A. Caraballo-Ortiz, “Mating system and fecundity of Goetzea elegans (Solanacea) an endangered tree of Puerto Rico,” in Department of BiologyUniversity of Puerto Rico-Rio Piedras, University of Puerto Rico, San Juan, Puerto Rico, 2007.

G. Allen-Wardell, P. Bernhardt, R. Bitner et al., “The potential consequences of pollinator declines on the conservation of biodiversity and stability of food crop yields,” Conservation Biology, vol. 12, no. 1, pp. 8–17, 1998.

A. Abebe, Market chain analysis of honey production: in atsbi wemberta district, eastern zone of tigray National Regional State [M.S. thesis], Department of Agricultural Economics, School of Graduate StudiesHaramaya University:, Alamaya, Ethiopia, 2009.