Mapping the Population Density of Managed Honey Bee (Apis Mellifera) Colonies in Ontario, Canada: 2018

The number of managed honey bee, Apis mellifera, colonies in Canada has more than doubled in the past century (Statistics Canada, 2018; Agriculture and Agri-Food Canada, 2019). As greater numbers of colonies populate the country, it is important to understand the implications of colony density, as well as map and monitor the distribution of said colonies to provide a bearing for researchers, government officials and beekeepers.

In the province of Ontario alone, approximately 100,000 colonies are managed by, on average, about 2,900 beekeepers (OMAFRA, 2019; Agriculture and Agri-Food Canada, 2019). A mix of commercial (≥50 colonies) and hobbyist (<50 colonies) beekeeping operations comprise this population. The province occupies a landmass of over one million square kilometers, but the majority of this land is inhospitable to humans, resulting in approximately 94% of the human population occupying only around 14% of the provincial landmass (Statistics Canada, 2021). Inherently, managed honey bees exist in correlation with human populations and therefore the 100,000 colonies in Ontario would also be expected to occupy a small fraction of the provincial landmass. However, population density distribution metrics of honey bee colonies are not as well documented as that of humans.

Beekeeping is a landless branch of animal husbandry and as such, there is limited control over the interactions of these managed bees and the surrounding environment. Therefore, the beekeeping community is heavily intertwined and reliant on one another to prevent disease spread and minimize the pool of bacteria and viruses in a given region. Therefore, knowledge about the geographic colony density distribution is of critical importance.

Host population density is a key component in basic host-parasite ecology (Kilpatrick & Altizer, 2010). In a density-dependent disease transmission model, the probability of a susceptible individual becoming infected is dependent on the number of other hosts in the same area that are already infected. This population density dependence is true for honey bees and the transmission of numerous viral and bacterial pathogens (Brosi et al., 2017).

Infectious disease modeling has demonstrated that the speed at which a disease spreads is a function of the contact rate (i.e., the rate of contacts sufficient to transmit a pathogen between infected and susceptible individuals) as well as agent (i.e., pathogen resilience and virulence), and environmental factors (i.e., climate and population density). Favourable environmental conditions, including the density of individuals around the susceptible host, and the homogeneity of individual mixing can accelerate the horizontal spread of disease (Lindahl & Grace, 2015). Increasing the density of hosts has the potential to increase the ratio of susceptible to infected individuals, thus increasing the contact rate and the likelihood of a successful transfer of the disease (Tarwater & Martin, 2001).

Humans have artificially increased the density of A. mellifera colonies as a means for industrializing agricultural farming and honey production. This is believed to have accelerated the disease transmission rates between colonies and thus increased the overall prevalence of various diseases and pests (Seeley & Smith, 2015). This increase in colony density, both within and between yards can be assumed to impact neighboring colonies and yards, with respect to their level of risk for disease and pests and influences their need for management strategies to mitigate their impacts.

Horizontal pathogen transmission between colonies of honey bees may occur by several mechanisms, namely, (i) robbing of honey stores, (ii) drifting, (iii) direct contact during foraging, (iv) contact with infected materials or food sources in the environment (Fries & Camazine, 2001), and (v) venereal diseases spread through mating (Chen & Siede, 2007). All of which may conceivably be influenced by the density of colonies in a region due to (i) greater competition during foraging, leading to intensified forage and robbing behaviour if nectar becomes limited, and (ii) higher inter-colony contact rates by other such means as drifting or fomites.

Although increased colony density may have significant implications on the behaviours of honey bees and the transmission of pests and diseases between colonies, an absence of information on local colony density values across the province of Ontario, Canada, prevents these associations from being investigated. Crude measures of density offer minimal practical value as they do not take into account heterogeneity or patterns of colony distribution. Geostatistical kriging methods offer a statistically sound means to convert regional counts of colonies (derived from registration data) into spatially continuous estimates of density for a larger study area (Matheron, 1963; Krige, 1981).

Detailed estimates of colony density provide a basis for epidemiological studies and practical disease intervention efforts, while routine colony inspection could also benefit from colony density estimates to implement a sampling strategy informed by the population distribution (i.e., assigning inspections to higher density locations). Ecological research may also in time reproduce such estimates to monitor the population growth and changes in distribution over multiple years. Moreover, beekeeper and community education efforts can benefit from graphic visualization of data patterns otherwise hidden in tabular presentations. Such graphics and maps regarding colony density can be part of the knowledge on mobilization tools, which are commonly used to bridge the gap between research and practical uses of results through the dissemination of research findings in an accessible form for broad audiences (Government of Canada, 2019).

While tabulated data is sufficient for basic data analysis, more advanced statistical mapping and modeling applications require data interpolated to grid locations across a study area. Estimating the spatially continuous colony density for the province allows for future research projects to extract data as raster or grid data at the necessary resolution. Because maps of the population density are of interest to both researchers and industry professionals, an objective of the present study is to develop an interactive map depicting regional colony density that serves as an online knowledge mobilization tool between beekeepers and government agencies.

The general goal of this study was to provide an accurate depiction of honey bee colony density in Ontario and identify spatial trends or patterns in the distribution of colonies. The specific objectives of this study were (1) to estimate the local density of managed A. mellifera colonies across Ontario through geostatistical interpolation of regional registry data, (2) to map the estimated density values as a density heatmap for applications requiring spatially continuous or high detail information (i.e. GIS analysis requiring raster data), (3) and to design an online knowledge mobilization tool of crude regional density values for applications requiring simple regional density metrics (i.e. public outreach or education requiring intuitive and interactive tools).

Records of registered bee colonies in Ontario were provided by the Ontario Ministry of Agriculture, Food and Rural Affairs (OMAFRA). The received records consisted of all registered bee yards in 2018, which was the most current data at the time of analysis. This data consists of the location to which a yard is registered by the apiarist and therefore assumes colony stationarity.

To maintain the privacy of beekeepers, all bee yard locations and their colony counts were aggregated to the Census Consolidated Subdivision (CCS) level before being released by OMAFRA. The CCS data divides the province of Ontario into 273 regions of varying sizes. The data consists of the total number of colonies and the total number of bee yards per CCS. A boundary file of the CCS subdivisions was obtained from Statistics Canada (Statistics Canada, 2019).

Within R (R Core Team, 2019), the boundary file was projected using the Lambert Equal Area (LEA) projection, because of the ability to retain the integrity of a square kilometer unit over a large geographical body, which is important for population density mapping. Bee colony densities per square kilometer were calculated and joined to the geographic centroid of each CCS in Ontario. The colony density values were heavily skewed and therefore a log-normal kriging approach was used.

Kriging is a geostatistical process that allows the production of a spatially continuous surface from distinct sampling points. In this case, a spatially continuous estimate of the population density can be estimated from a collection of regionally aggregated colony registration counts.

The process of geostatistical kriging requires that a variogram model be fit to the data to determine the overall spatial dependence and the degree to which that dependence diminishes over a given distance. A spherical variogram model was fit to the data through maximum likelihood estimation (Ribeiro & Diggle, 2018). The range of the variogram was restricted to a maximum of 500 km, assuming local but not global stationarity. The model fit was assessed via the Akaike information criterion, as well as a visual inspection of a normal QQ-plot of the model residuals.

Ordinary log-normal kriging interpolation was performed, assuming a constant unknown mean across the study area. The process of log-normal kriging includes an automated procedure to transform the data and back-transform the prediction values during kriging. The R package, geoR was used for all spatial data management and analysis (Ribeiro & Diggle, 2018).

Isopleth mapping based on the back-transformed results from geostatistical kriging was used to illustrate the provincial population density of colonies with the use of the CCS data (Carrat & Valleron, 1992; Berke, 2004). All spatial regression model parameters including the variogram were estimated through maximum likelihood estimation.

Ordinary log-normal kriging interpolation was performed at the provincial level first and then separately applied to the Southern Ontario region following the removal of data from Northern Ontario. This was deemed appropriate due to Southern Ontario's distinctively smaller regions, which are difficult to recognize on a provincial map. Furthermore, Northern and Southern Ontario differ dramatically in terms of size and crude colony density, so it cannot be assumed that a single spatial dependence structure accurately represents the entire province. The two population density maps for Ontario and Southern Ontario differ in their colour scales and are to be compared with this in mind.

Recognizing that there is a shared need for colony density information by researchers and industry professionals, a practical knowledge mobilization dashboard was developed. Using the leaflet package in R (Cheng et al., 2018), an interactive choropleth map at the CCS level was produced. The map allows for the exploration of the colony densities across the CCSs with the ability to zoom to finer resolutions. Upon selection of a CCS region, a pop-up window with further colony density information is displayed. This internet dashboard can be updated upon the availability of new data.

The Ontario provincial dataset represents colony counts from 273 census consolidated subdivisions. The total number of colonies in Ontario in 2018 was n=87,777. A map representing the 273 CCSs and the associated crude density values is provided in Fig. 1.

Fig. 1

The 273 CCSs ranged in areal size from about 52km² to 446,000 km² [mean: 3,613±2,930 km²]. The crude colony density values for each CCS varied from 0 to 14.7 colonies/km², with an average crude colony density of 0.9 colonies/km² [0.88±1.75]. The colony density distribution is visualized with a choropleth map in Fig. 1.

The estimation of the semivariogram, using the spherical model is presented visually in Fig. 2, juxtaposed with the empirical semivariogram of the colony density data. The spherical variogram model was determined to best represent the spatial dependence structure of the data based on visual analysis of the variogram and the shape of the spherical function, as well as the AIC when compared to other models (i.e., exponential, gaussian, or linear). The range of spatial dependence based on the spherical model fit to the data is estimated at 360 km (Tab. 1). The estimated variogram shows a reasonably good fit between the model and the empirical variogram (Fig. 2a), and the standardized model residuals moderately maintain the assumption of normality with deviations at the tails (Fig. 2b). An isopleth map of colony density for Ontario, derived from geostatistical kriging of the crude regional density values, is shown in Fig. 3. The average colony density across the province was estimated at a level of 0.373 colonies/km². Location-specific colony density estimates ranged from 0 to 8.93 colonies/km². An uneven distribution of colonies appears evident based on the visual examination of the estimated density map. The colony density map provides evidence of a heterogeneous distribution of colonies in Ontario, with a majority of managed honey bee colonies residing in the south of the province (<−500N, see Fig. 3), and more specifically in Niagara and surrounding regions (500E, −900N, see Fig. 3).

Fig. 2

Fig. 3

Similar to the provincial data analysis, the colony densities for Southern Ontario CCS regions were positively skewed and thus log-transformed for variogram model fitting. This transformation was found to be successful in improving the fit of the empirical variogram estimation based upon visual inspection of a histogram. Like the previous variogram estimation, a spherical model was found to be appropriate for kriging purposes (Fig. 4a). The model residuals appear to better maintain the assumption of normality, compared to the complete province model, with an approximate mean nearing 0, based on visual inspection of the quantile-quantile plot (Fig. 4b). A practical spatial dependence range of 191 km was estimated for Southern Ontario by the variogram model (Tab. 1). Kriging-estimated colony density estimates for Southern Ontario show a large proportion of colonies residing in the Niagara Peninsula (500E, −900N) with several less concentrated pockets in the Bruce/Grey (360E, −750N), Simcoe (450E, −760N), Wellington (425E, −850N) and Northumberland regions (600E, −760N), relative to other regions of Southern Ontario (Fig. 5). The average colony density for Southern Ontario was estimated at 0.863 colonies/km², with grid estimates ranging from −0.023 (0) to 10.4 colonies/km². This contrasts with the average colony density for Northern Ontario independently at 0.079 colonies/km².

Fig. 4

Fig. 5

Choropleth and isopleth mapping of the colony densities at a regional level suggests a disproportionate distribution of colonies with higher densities in the Niagara region and lower densities in the northern regions.

A screenshot of the interactive choropleth map is presented in Fig. 6.

Fig. 6

The results of this study illustrate the spatial distribution and density of managed honey bee colonies in Ontario. The mapped colony density estimates indicate an overall heterogeneous distribution across the province. When compared to the province as a whole, Southern Ontario displays a high density of colonies with an average density of around 8.6 colonies/10km². compared to Northern Ontario with a density of 0.8 colonies/10km², and the collective provincial average of approximately 3.7 colonies/10km². Because this study was focused on the managed honey bee population, this distribution of colonies is not surprising, as the Southern Ontario region comprises nearly 94% of the province's human population (Statistics Canada, 2017). However, this result may be surprising in the sense that Southern Ontario only accounts for 32% of the province's agricultural farm acreage (OMAFRA, 2020) where a sizeable proportion of commercial honey bees would be expected.

The large variation in colony density across the province gave rise to a substantial positive skewing of the data. This positive skewing is common in spatially dependent observations and has been shown to cause a proportional effect in the variogram estimates, where variation is dependent on the magnitude of the value (Manchuk et al., 2009). In this study, despite the skewing of data and the differing sill values, a proportional effect was not detected between the Ontario and Southern Ontario variogram models as the range was also found to differ. However, the skewness of the data is thought to be the primary reason for the lack of fit of the provincial variogram model.

Neither the variogram model estimation nor the residual QQ-plot for the whole province (Fig. 2) showed an exceptional fit to the data, and therefore the decision was made for Southern Ontario to be analyzed independently. By area, Northern Ontario comprises approximately 86% of the province but managed honey bee colonies in the region are sparse. The stark contrast between the colony densities in the two regions justified the decision for individual analyses since modeling the whole province with a single spatial dependence structure, while possible, is not ideal for accurately estimating model parameters. However, it was not possible to produce an independent model for Northern Ontario because the sample size of CCSs is too small without borrowing information from Southern Ontario. Consequently, the province-wide analysis was maintained to offer some representation to northern beekeepers.

When Southern Ontario was analyzed separately, there was less variation and skewness in the regional colony density values, and the estimated variogram was found to represent the data more closely (Fig. 4a). This improved model fit is further evident by the QQ-plot of the standardized residuals (Fig. 4b).

The Niagara Peninsula is most densely populated with bee colonies, with over five colonies/km² at minimum, and as high as over ten colonies/km² in the areas adjacent to New York State, USA (Fig. 5). Several factors may be at play for the magnitude of apiculture activity in Niagara. First, and most evident, is the prominent agricultural landscape of Niagara, and the need for commercial pollination to aid efficient crop production. Many of Niagara's top production crops are either self-pollinated, like peaches, cherries, plums, and grapes, or wind-pollinated such as grains, and therefore not insect pollinator-dependent. However, some of these self-pollinating crops are still commercially pollinated by honeybees to increase the yield and quality of the fruit (Klatt et al., 2014). Therefore, the abundance of agriculture may offer a limited explanation for the high density of honey bee colonies in this region. A second factor for the high bee colony population density in Niagara is the presence of commercial beekeeping operations and registered addresses for numerous commercial pollinator operations servicing farmland outside of Niagara or Ontario.

Niagara, being of the most southern regions in Canada, is a preferred home for commercial beekeepers, in part due to the relatively mild winters offering favourable conditions for overwintering and an early spring harvest, allowing for colonies to build up strength in time for eastern Canada blueberry pollination in May and June.

Some spatial misclassification may result from the nature of the government data used for this study. The colony location data indicate the locations where colonies are registered, typically their overwintering site, but may not always be accurate to the locations where the colonies reside during the entire season. These registrations and resulting data represent a snapshot in time and the numbers of colonies and their locations can change throughout the season. This issue becomes more dramatic when out-of-province pollination demands, such as blueberry pollination in eastern provinces, significantly reduce the number of colonies residing in Ontario, as roughly 20–35% of colonies leave at this time for pollination in New Brunswick, Prince Edward Island, and Quebec (OMAFRA, 2013). This considerable movement of honey bee colonies causes the population density to be a dynamic and complex topic of study.

Another study limitation is the potential for underreporting of colony ownership. Under-reporting can occur when beekeepers do not accurately declare ownership of colonies to the provincial government (OMAFRA) and instead report fewer than the actual number of colonies, or when beekeepers do not report keeping bees at all. Assuming that beekeepers have no regionally varying motivation to underreport, the underreporting effect should be non-differential and affect all areas proportionally. Thus, the spatial pattern would not be biased, but the colony density estimates have the potential for minimal underreporting. Under the Ontario Bees Act, “No person shall be a beekeeper in Ontario without a certificate of registration issued by the Provincial Apiarist” and this registration includes the disclosure of the number of colonies in possession (Ontario Bees Act, 1990). Thus, to the knowledge of the Ontario Ministry of Agriculture, Food and Rural Affairs, along with the provincial beekeeper, the data used in this study are accurate and well representative of the population of managed colonies in Ontario.

To protect the privacy of beekeepers and yard locations, this study only received access to regional data at the level of the CCS. The use of aggregated data inherently results in some loss of information for data analysis because the exact spatial distribution within each region is averaged. Regions with larger geographic areas are more likely to be affected by this regional aggregation. This is expected to be a larger issue in Northern Ontario where regions with a large landmass accounted only for a small number of colonies. For Southern Ontario, the bias from the regional aggregation of the data is considered relatively small because colonies may not be tied to a single point but moved within the region during the season. Therefore, point data would not necessarily be more accurate in regions that are small relative to the movement of colonies. Yard locations may be updated by an inspector at the time of inspection if discrepancies are found.

If this study could be based on precise locations of the colonies (i.e. spatial point data) the findings could change, which is known as the modifiable areal unit problem in geographic epidemiology (Waller & Gotway, 2004), as clustering within the CCS is likely to exist in yard-level and operation-level spatial groupings of colonies. This would be expected to impact the more remote regions of Northern Ontario, where the few existing colonies are likely operated by only a handful of beekeepers in a limited number of yards. This clustering hierarchy of yards, operations, and CCS regions contributes further complexity to the topic of honey bee population density when reliant on aggregated data. In Southern Ontario, where regions are smaller and the number of colonies is greater, the variogram model accounts for spatial correlation of colony locations and kriging can adequately accommodate and exploit these spatial correlation structures to produce continuous spatial estimates. The clustering of colonies is a spatial data pattern related to the location of placed colonies by beekeepers and does not necessarily relate to the biological flight range of bees.

From the results of this study, it follows that the density of managed honey bee colonies in select areas of Ontario is quite high (upwards of 10 colonies/km² in areas of Southern Ontario), especially when compared to the apparent natural density of feral honey bees in forests of New York at roughly 1 colony/km² (Seeley, 2007; Seeley et al., 2015). While the results from this study draw no conclusions on the health status of bees in high-density regions of Ontario, they do identify areas that may be of interest for future studies based on the known implications of population density on pathogen and pest transmission between colonies. Future work on the spatial distribution of diseases among managed honey bee colonies could consider the findings from this study to standardize disease counts based on the population-at-risk to prevent misrepresentation of disease frequency in areas of high population density.

Based on this study, Niagara is hypothesized to possess a greater total prevalence of communicable diseases due to the higher number of colonies at risk, but disease prevalence is also dependent on management practices and treatment schedules. In instances of Varroa destructor-induced colony death, an abundance of colonies in close proximity to the weakened colony may increase the probability of robbing to occur and may transfer mites between colonies (Frey & Rosenkranz, 2014; Peck & Seeley, 2019). Furthermore, the transmission of Paenibacillus larvae by means of robbing may occur more readily in the higher density regions identified, where intercolonial distance is within the 1 km range, which is a reported risk factor (Lindström et al., 2008). These high-density regions may also selectively support more highly virulent strains of P. larvae that can induce more severe clinical symptoms without concern of exhausting available hosts, where pathogens in low-density regions must be more sparing of their hosts and not be lethal enough to out-pace host reproduction (Fries & Camazine, 2001). Horizontal transmission by means of drifting could also be expected to be elevated in regions of high density such as the Niagara peninsula based on previous literature on the consequences of density on drifting (Seeley & Smith, 2015). Furthermore, a region possessing a high density may not be indicative of close proximity between neighbouring yards as the number of colonies per yard can vary.

The results of this study provide insight into the distribution of colony density and distribution across Ontario, Canada. The overall population density structure across the province is heterogeneous with two main patterns. Southern Ontario is highly concentrated compared to Northern Ontario, and the Niagara Peninsula is highly concentrated compared to the remainder of Southern Ontario. Future studies should aim to ground-truth these results, such as using prevalence data of honey bee pathogens to confirm an association between high colony density and high prevalence of various pathogenic diseases and pests (e.g., varroa mites). Results from this study can inform risk-based apiary health inspections through a focus on areas at potentially greater risk to transmit disease. Furthermore, this study highlights the methodologies available for similar mapping projects to be conducted elsewhere, including areas with non-uniform densities. Disease transmission is complex and multifaceted, and therefore, it should not be concluded from this study alone that the areas identified as high density possess higher than average disease or pest risk. These findings are meant to expose the potential implications of high colony density and give beekeepers, researchers and government officials a statistically sound estimation of the densities in all regions of Ontario. Furthermore, visually intuitive and interactive maps developed as part of this project provide adequate information distribution to members of the commercial and hobbyist beekeeping community.