Agricultural Landscape and Pesticide Effects on Honey Bee (Hymenoptera: Apidae) Biological Traits.

Sixteen honey bee (Apis mellifera L.) colonies were placed in four different agricultural landscapes to study the effects of agricultural landscape and exposure to pesticides on honey bee health. Colonies were located in three different agricultural areas with varying levels of agricultural intensity (AG areas) and one nonagricultural area (NAG area). Colonies were monitored for their performance and productivity for one year by measuring colony weight changes, brood production, and colony thermoregulation. Palynological and chemical analyses were conducted on the trapped pollen collected from each colony and location. Our results indicate that the landscape's composition significantly affected honey bee colony performance and development. Colony weight and brood production were significantly greater in AG areas compared to the NAG area. Better colony thermoregulation in AG areas' colonies was also observed. The quantities of pesticides measured in the trapped pollen were relatively low compared to their acute toxicity. Unexplained queen and colony losses were recorded in the AG areas, while colony losses because of starvation were observed in the NAG area. Our results indicate that landscape with high urban activity enhances honey bee brood production, with no significant effects on colony weight gain. Our study indicates that agricultural crops provide a valuable resource for honey bee colonies, but there is a trade-off with an increased risk of exposure to pesticides.

Pesticides are thought to be a principal factor causing honey bee decline in addition to the parasitic varroa mite (Varroa destructor; Johnson et al. 2010, Van der Sluijs et al. 2013). The honey bee (Apis mellifera L.) is one of the most efficient pollinators of a wide range of plants and crops (Robinson et al. 1989, Morse and Calderone 2000, Calderone 2012, Klatt et al. 2014). Most North American crops, as well as crops worldwide, are able to produce seeds and fruits without pollination. However, bee pollination was proven to increase the value of some fruits such as strawberries through higher yields and better quality (Klein et al. 2007, Klatt et al. 2014). The number of honey bee colonies in the United States has nevertheless declined by 45% over the past 60 yr (NAS 2007), not just because of agrochemical exposure (Johnson et al. 2010, Van der Sluijs et al. 2013) but also a result of various pathogens, parasites (Fries 2010, Dainat et al. 2012), and other factors such as changing farm demographics. Moreover, some environmental factors play a central role in colony losses, such as habitat loss or changes, poor nutrition, inadequate foraging flora, and the transportation stress induced by the excessive transhumance of honey bee colonies to provide pollination services (Naug 2009, Ahn et al. 2012, Simone-Finstrom et al. 2016).

Pesticides are necessary for pest control in agriculture, increasing crop production, and providing worldwide food security (Potts et al. 2010). However, many insecticides are toxic to bees at very low doses (Williamson et al. 2014, Chakrabarti et al. 2015), and they may cause significant disorders at sublethal doses in colony dynamics and the division of labor of honey bee colonies (Mackenzie and Winston 1989) by affecting honey bee behavior, orientation, dance communication, and return flights (Vandame et al. 1995, Fischer et al. 2014, Williamson et al. 2014). One other negative aspect of the pesticides, more precisely the neonicotinoids, is their capability of supressing honey bee immunocompetence that might lead to an impaired disease resistance capacity (Di Prisco et al. 2013, Brandt et al. 2016). In other words, bee-toxic pesticides can threaten bee health, particularly if misused or applied without consideration to effects on honey bees or pollinators.

Honey bees are not only facing threats of agricultural pesticides, but also those used in-hive by beekeepers themselves to treat pathogens and parasites. For instance, beekeeper-applied pesticide residues at various concentrations were identified in honey bee products and hives (Lambert et al. 2013), including a nearly ubiquitous detection of coumaphos and fluvalinate that are used to control varroa mites (Mullin et al. 2010, Vanbergen and The Insect Pollinators Initiative 2013, Ravoet et al. 2015). The long-term accumulation of these chemical residues, their combination, and their synergistic effects with other health factors (e.g., pathogens) could eventually lead to unforeseen negative effects on colony survival (Chauzat et al. 2006, Lambert et al. 2013).

The performance of a honey bee colony is the result of complex dynamics involving many internal and external factors. In order to thrive or even just survive, a colony must withstand various environmental conditions, including parasites and pathogens, and reach an optimal balance in its food supply or reserve and population size (Adam 1983, Ruttner 1988, Buchler et al. 2014, Meixner et al. 2014). Colony decline or loss may sometimes be linked to one main factor, but is more often the result of a gradual deterioration in the population health and immunity because of multifactorial synergistic effects (Alaux et al. 2010, Nazzi et al. 2012, Di Prisco et al. 2013).

Two biological traits (colony weight and brood production) are the main characteristics that can fairly assess colony performance in Terms of honey production (nectar collection and environmental flora) and population size (queen health and productivity; Louveaux 1973, Alburaki et al. 2015, Meikle et al. 2016). Additionally, other important factors are also involved in colony performance, such as genetics and pathogenicity (Tarpy et al. 2013). Thermoregulation of the honey bee nest cavity is another important factor for colony survival, and it is regulated by the worker bees with high precision (Winston 1987). It has been determined that even a small deviation from the optimal brood-nest temperature of 35 °C (Seeley and Heinrich 1981, Southwick 1985, Jones et al. 2004) can significantly influence brood development and the resulting production of adult bees (Matthias et al. 2009). Moreover, adult bees reared at suboptimal temperatures are more vulnerable to certain pesticides (Medrzycki et al. 2009).

In this study, we measured three key elements of honey bee colony health—colony weight, brood production, and colony thermoregulation—in different landscapes and with different risks of pesticide exposure. We then evaluated honey bee colony performance in replicated exposure groups in an effort to tease apart the relative effects of pesticides and environment on colony health.

Materials and Methods

Honey Bee Colonies

This experiment was conducted on 16 honey bee colonies. In May 2015, new colony divisions—each equivalent in size—were made from existing, overwintered colonies of standard stock. Sixteen Carniolan (Apis mellifera carnica) queens, artificially inseminated, were purchased from a commercial queen producer at the same period and introduced into the new colonies. The new queens and divisions were housed in new 10-frame Langstroth hives (Table 1).

Table 1

Background, geographical coordinates, and distribution of the experimental apiaries. GIS landscape classification of each apiary location is calculated on a 2.5-km radius. Main crops and pesticides used in the studied locations as well as beekeeping procedure applied on experimental hives during this study

	Apiary 1	Apiary 2	Apiary 3	Apiary 4
Location	Jackson	Milan	Yum-Yum	Chickasaw
Code in Graphics	J	M	Y	C
GPS Coordinates	35° 37′38.81″ N	35° 56′21.03″ N	35° 21′17.00″ N	35° 23′15.63″ N
	88° 51′03.55″ W	88° 43′12.64″ W	89° 20′50.73″ W	88° 46′56.98″ W
Landscape classification	Low AG area with urban activity	Moderate AG area	High AG area	NAG area
Main crops	Corn, Soybean, Sorghum, Cotton			None
Pesticide used	Neonicotinoids, organophosphates, pyrethroids, and other insecticides; fungicides, herbicides			None
Pesticide application	Coated-seed and foliar			None
Period of study	May-2015 to April-2016
No. of colony	4	4	4	4
Colony code	H1, H2, H3, H4	H5, H6, H7, H8	H9, H10, H11, H12	H13, H14, H15, H16
Type of hive	10-frame Langstroth
Queen origin	2015—Apis mellifera carnica
Hive equipment	Temperature sensor (fixed on frame N°5)
	None—Manually weighed		Electronic scales
Weight data	May 2015 to March 2016
Brood data	May to September 2015
Nest temperature	May 2015 to March 2016
Trapped pollen collected	8 times/hive, June to September 2015
Chemical analysis/trapped pollen	4 samples	4 samples	4 samples	4 samples
Super added	One/H4	One/H8	One/H9 and H10	One/H13
Winter feeding	Four times/1:1 Sugar syrup only
Total colony perished	2	1	0	2

Location and Landscape

The four experimental locations (Jackson, Milan, Yum-Yum, and Chickasaw) were carefully chosen based on potential exposure to pesticides and cropping intensity (Fig. 1). As organic crop production is rare in western Tennessee, the control treatment location (Chickasaw) was placed in a NAG area in a Tennessee state park where no cropping areas were in the near vicinity (Fig. 1). Geographical information system (GIS) studies were conducted on a 2.5 km-radius from each candidate location. The GIS analysis was conducted in order to assess the AG areas or crop fields available for honey bees within a typical foraging distance for honey bees (Seeley 2010). The GIS study was performed using Esri ArcGIS software (Redlands, CA).

Fig. 1

Geographical location of the four studied apiaries in western Tennessee, USA. The GIS study to determine the landscape nature at each location was conducted on a 2.5-km radius. The landscape classifications included Jackson (low AG area with urban activity, 19%), Milan (moderate AG area, 55%), Yum-Yum (high AG area, 71%), and Chickasaw (NAG area, 5%).

Weight Evaluation

Two different methods were used to weigh the experimental hives. First, hives were weighed manually by using a scale with a 0.01 kg-sensitivity. Second, hives were weighed using electronic scales on which hives were permanently kept until the end of the experiment in March 2016. In the first case, colonies of two locations (Jackson and Milan) were weighed biweekly from May 2015 to March 2016 (Table 1). The electronic scales used at Yum-Yum and Chickasaw provided near-continuous measurements of hive weights. Data related to the set up and use of the electronic scales are available in a recently published study (Meikle et al. 2016).

Brood Production

Brood development was evaluated biweekly for each colony from May to September 2015 (Table 1). Frames containing capped brood worker were photographed to quantify the worker production of each colony. Frames were temporarily removed from a hive, brushed to remove adult bees, placed onto a special apparatus equipped with a digital camera, and photographed on both sides. Photos were subsequently downloaded on computer, and worker capped brood cells were manually counted for each studied colony using IMAGEJ software (Abràmoff et al. 2004). Although time consuming, the brood counting method used in our study is highly accurate and provided near exact counts of capped brood in hives (Bertrand et al. 2015).

Inner Colony Temperature

Each of the 16 hives was equipped with a temperature sensor (iButton) fixed on the upper side of the middle frame of the hive. In order to avoid bees’ direct contact with the device, each sensor was covered with a 6- by 6-cm metallic screen (see DOI in supporting information). All iButton devices were set to perform one temperature read every hour with a sensitivity of 0.06 °C. These data were collected throughout the experiment going from May 2015 to April 2016 (Table 1) and were periodically collected by downloading them from the iButton directly to a computer.

Palynological Analysis

Trapped pollen was collected from each hive at eight different time points from June to September 2015 (Table 1). Pollen collected from each colony was desiccated at 37 °C for 48 h and stored individually at −20 °C. A detailed palynological analysis was performed on the trapped pollen to define its origin and composition. Briefly, each dried sample was homogenized and 1 g of each sample was used to determine the botanical origin of the pollen loads with ∼300 observed pollen grains per sample. The taxonomic diversity of pollen samples for each colony at each sampled date was determined by observing the total surface of slides (Loublier et al. 2003).

Pollen Pesticide Residues

Pesticide residues were quantified in the trapped pollen using liquid chromatography-mass spectrometry (LC-MS; Barnett et al. 2007, Walorczyk and Gnusowski 2009). All chemical analyses for pesticide residue detection were processed at the USDA National Scientific Laboratories in Gastonia, North Carolina. Trapped pollen of each month and apiary were pooled and a comprehensive chemical analysis was run for 16 pollen samples that included 174 chemical compounds or molecules.

Statistical Analysis

Colony weights of our four experimental apiaries were taken in kg except for hives on electronic scales (Yum-Yum and Chickasaw), for which Hobo devices provided the weight data in milliampere (mA). Data obtained in (mA) from Hobo devices were translated to kg as described in (Meikle et al. 2016). Data weights across this study were given as a net value of the bee population weight after omitting the weight of the hive box and any additional equipment added through the season such as supers and feeders. However, the four sugar feeds were considered a part of the bee population development. Concerning the brood data, capped worker brood counted from frames of each colony were summed for each colony and sampling date (Alburaki et al. 2015, Alburaki et al. 2016). Outlier values were not omitted from the datasets and were treated as such to fairly assess the studied factors. However, disturbances in the weight and colony temperature values resulting from hive management were discarded from the dataset before statistical analysis.

Statistical analysis and figure generation were carried out and generated in the R environment (R Core Team 2011). Variables of this study included—1) colony weight (kg), 2) brood production, 3) inner hive temperature (C°), and 4) apiary location. Data were statistically treated per location (four groups, each includes four biological replicates) to study the landscape and the putative impacts of the pesticide exposure on honey bee colonies.

One-Way ANOVA Tests

Variables were first tested for normality using the Shapiro-Wilk test and none of them was normally distributed (see DOI/Variable distribution—Shapiro test). We attempted to normalize our data by log-transformation, ran again Shapiro-Wilk tests and carried out Q-Q plots to visualize the data distribution. Log-transformation failed to normalize our dataset; nonetheless, analysis of variance (ANOVA) was carried out at a 95% confidence level. ANOVA is not very sensitive to moderate deviations from normality (Glass et al. 1972, Harwell et al. 1992, Lix et al. 1996), and our compared variables were equal in size. Correlations between colonies for each variable were performed using the R libraries “Performance Analytics” and “Corrplot”, respectively. Pearson method was used in all of the correlation analyses performed in this study.

Generalized Linear Mixed Analyses

In order to assess the accumulative effects of all variables together on honey bee colony health, linear models (lme) were used. Our data were not normally distributed and includes repeated measures; therefore, a general model was indicated. In order to fairly assess the response variables by considering the experimental uncontrollable factors (such as environment and flora variability), random effects must be considered in our models. Thus, the most appropriate model for our experiment design and the nature of our data was the “generalized linear mixed-effects model” (GLMM) in which we considered “fixed effects” as explanatory variables and “random factors” as uncontrollable experimental factors (Baayen et al. 2008). To improve accuracy in assessing the effect of the landscape and pesticide exposure on colony health without violating the independence assumption and accounting interdependencies, the factor “AG area” (1| AG) was treated as a random factor in all our models (Bolker et al. 2009). As our models included fewer that three random effects, and no binary response variables were tested, penalized quasi-likelihood (PQL) was used to approximate the likelihood to estimate the GLMM parameters (Schall 1991, Breslow and Clayton 1993, Wolfinger and O'Connel 1993) as explained in the model below. The GLM models used in this study required both “nlme” and “MASS’ packages (Bates et al. 2014) and can be summarized as follows:

Results

Landscape Classification

The locations and their agricultural classification based on the GIS were as follows (Fig. 1 and Table 1): Jackson (low AG area with urban activity), Milan (moderate AG area), Yum-Yum (high AG area), and Chickasaw (a natural park that contains essentially no agricultural activity; NAG area). Within a 2.5-km-radius foraging distance, honey bees had access to 19, 55, 71, and 5% of total agricultural landscape in the Jackson, Milan, Yum-Yum, and Chickasaw sites, respectively (Fig. 1). The remaining landscapes consisted of forest, woodland, open water, and urban habitat (e.g., buildings and roads). Chickasaw was considered the control treatment of this study with only 5% AG area. Jackson location was the only location with a relatively high urban activity (46%, Fig. 1).

Crop Fields

Based on GIS analyses and observation of the crop fields surrounding each apiary, four main summer crops were cultivated in the studied areas: 1) Corn (Zea mays), 2) Cotton (Gossypium hirsutum), 3) Soybean (Glycine max), and 4) Sorghum (Sorghum bicolor), which are typical for those grown in West Tennessee. Variable types of clover (Trifolium sp.) were also observed in all locations. Winter wheat (Triticum aestivum) was also commonly present in AG areas.

Colony Weight Development

Depending on the location, colony weights were periodically or continuously monitored from May 2015 to March 2016 (Fig. 2). Disturbance in the weight data of the electronic scales because of hive examination, sampling, or measurement was discarded from the dataset before conducting the statistical analyses. Overall, colony weight averages of each apiary were: 1) Jackson = 6.99 ± 0.37 kg, 2) Milan = 10.28 ± 1.04 kg, 3) Yum-Yum = 12.29 ± 0.03 kg, and 4) Chickasaw = 6.49 ± 0.01 kg (Fig. 3A). Milan and Yum-Yum weights were significantly greater than Chickasaw (P < 0.001), while Jackson was not different than Chickasaw (Fig. 4). The weight correlation matrix revealed three different colony scatters, 1—colonies negatively correlated (bottom left in red), 2—colonies positively correlated (bottom right in blue), and 3—colonies in between (middle branch), Fig. 3A.

Fig. 2

Colony weight records taken manually from May 2015 to March 2016 showing the net weight (kg) of hives located in Jackson (A) and Milan (B). ANOVA revealed significant differences in weight development for H8 only (P < 0.001).

Fig. 3

(a) Weight development, (B) brood production, and (C) inner colony temperature exposed by location (treatment). Location codes are C = Chickasaw, J = Jackson, M = Milan, and Y = Yum-Yum. Outlier values were kept in the dataset for higher accuracy, significant differences among groups are indicated by *** (P < 0.001) and ** (P < 0.01). Heat maps are the correlation matrixes of overall colony weight, brood, and temperature.

Fig. 4

Brood production for each studied colony in 2015. Data are log-scaled and exposed by date and colony, data represents the total number of sealed worker brood in the hive as determined by manual counting using IMAGE J software.

Among the eight sampling dates, the highest brood production was recorded in Jackson’s apiary (12,840 capped workers on the month of July) followed by Milan’s apiary (11,996) in September (see DOI in supporting information). Chickasaw’s colonies were consistently lower than other locations in brood production and did not exceed 8,031 capped worker cells at any time point. Disruptions and fluctuations in brood rearing were more obvious in Chickasaw than other locations (Fig. 4). Brood production average per location were as follows: 1) Jackson = 8,987 ± 820, 2) Milan = 8,348 ± 600, 3) Yum-Yum = 7,486 ± 707, and 4) Chickasaw = 5,568 ± 481 worker capped cells. Our brood data indicate two queen losses in our experimental colonies, both from the AG areas including H2 on 24 August and the H9 on 14 September (Fig. 4). A late start in brood production was easily noticeable for H3, which did not produce brood until 15 June (Fig. 4), likely because of another requeening event.

ANOVA revealed significantly higher brood production in Jackson and Milan (AG areas) than in Chickasaw (NAG area), (P < 0.01; Fig. 3B). The correlation matrix clearly revealed two correlation types, negative and positive (Fig. 3B). Most colonies tended to positively correlate with each other in their brood development (blue) with exception of H2 and H9 that showed negative correlations (red) with all other experimental colonies (Fig. 3B). H15 and H16 (from the NAG site) were the two most neutral colonies for brood production (Fig. 3B).

Inner Beehive Temperature

In total, iButton devices provided 10,224 temperature-reads per colony during our one-year study except for colonies that perished. Periodic disturbances in the temperature related to our sampling efforts were omitted from the dataset. Colonies showed better ability to adjust their 35 °C colony temperature from June to October 2015. Greater oscillations were noticeable during winter (from November 2015 to late February 2016; Fig. 5). Overall averages of inner hive temperature for colonies at each location were as follows: 1) Jackson = 31.75 ± 0.03 °C, 2) Milan = 31.22 ± 0.03 °C, 3) Yum-Yum = 31.04 ± 0.03 °C, and 4) Chickasaw = 29.69 ± 0.03 °C. AG areas’ colonies maintained an overall colony temperature of 31.33 ± 0.03 °C versus 29.69 ± 0.03 °C for colonies in the NAG area. ANOVA indicated significantly higher hive temperature in AG areas (Jackson, Milan, and Yum-Yum) versus NAG area (Chickasaw) (F= 1,033, n = 154,776, P < 0.001; Fig. 3C). Concerning thermoregulation, the correlation matrix mostly showed significant positive correlations (blue) among colonies (Fig. 3C) with some exceptions, such as H10 vs. H2 (r = −0.87, P < 0.001) and H9 vs. H2 (r = −0.69, P < 0.01).

Fig. 5

One-year (2015–2016) monitoring of internal beehive temperature (°C) for each of the 16 studied colonies. Temperature was measured using iButton device fixed on the middle frame of each hive. Punctual perturbations because of hive examination and sampling were omitted from the dataset.

Colony Mortality

In total, five of the 16 colonies perished during this one-year study, with two colonies located in Jackson (H2, H3), two in Chickasaw (H13, H16), one in Milan (H6), and none in Yum-Yum (Table 1). In Jackson, colony H2 perished in November 2015 followed by H3 in February 2016. Visual examinations indicated that queen loss and failure to requeen was the main reason why H2 died. H3, however, likely died from typical winter mortality (i.e., starvation). In Milan, colony H6 died in February 2016 with a small bee population size that was unable to thermoregulate the hive’s cavity (Fig. 5). The most obvious and remarkable colony losses were those of H13 and H16, both of which perished simultaneously at the end of our experiment (March 2016) likely because of starvation (see DOI in supporting information). Curiously, those colonies had the highest weight and brood production at that location across the year (Fig. 5, 6B).

Fig. 6

Outputs of the eight electronic scales showing colony weight development from May 2015 to March 2016 of two locations: (A) Yum-Yum and (B) Chickasaw. Missing data in (A) are because of battery malfunction of the Hobo device. ANOVA revealed significant differences in weight of both hives H9 and H13 (P < 0.001).

Pollen Identification

Various types of pollen were microscopically identified in the trapped pollen of each hive (Table 2). All pollen grains of the main four crops cultivated around the experimental colonies were identified in most samples with the notable exception of cotton. Cotton pollen was not identified in any hives (Table 2). Corn pollen, normally considered as not attractive to honey bees, was only identified at substantial levels at the Jackson (12.6%) and Yum-Yum locations (12.3%). Soybean pollen was commonly collected at the Milan (46.3%) and Yum-Yum (32.9%) apiaries. Sorghum pollen was also commonly collected at Milan (12%) and Yum-Yum (14.1%). In the NAG area (Chickasaw), no soybean pollen was identified, and less than 1% of the pollen identified was either sorghum or corn (Table 2).

Table 2

Palynological analysis of the trapped pollen collected by forager bees for each studied hive and location. Percentage of the crop pollen identified in each hive as well as the four major noncrop pollen. Pollen were identified under microscopy and percentage calculated based on ∼300 counted grains per slide per sample (Loublier et al. 2003)

Apiary location	Colony	Corn	Soybean	Sorghum	Cotton	Clover Trifolium sp.	Elderberry Sambucus spp.	Sumac Rhus typhina	Plantain Plantago sp.
Apiary 1 (Jackson, Low AG)	H1	17.6	2.7	4.5	–	5.2	20.2	27.7	16.5
	H2	9.1	2.1	–	–	21.4	1.4	24.2	10.6
	H3	20.7	0.8	0.3	–	7.0	13	28.1	31.5
	H4	3.0	1.4	4.7	–	16.2	3.4	16.6	11.1
Avg		12.6	1.7	2.3	–	12.4	9.5	24.1	17.4
Apiary 2 (Milan, Moderate AG)	H5	0.7	37.6	17.1	–	3.7	–	57.5	2.0
	H6	0.3	48.2	–	–	6.6	41.6	0.6	2.6
	H7	0.3	52.5	–	–	13.7	–	0.6	15.9
	H8	3.0	47.0	30.7	–	22.5	0.3	3.2	22.2
Avg		1	46.3	11.9	–	11.6	10.4	15.5	10.7
Apiary 3 (Yum–Yum, High AG)	H9	–	23.1	3.4	–	48.4	50	32.6	0.5
	H10	–	33.2	11.7	–	24.7	–	64.3	0.3
	H11	–	42.4	21.7	–	0.3	–	–	–
	H12	49	33	19.8	–	32.4	–	16.2	–
Avg		12.2	32.9	14.1	–	26.4	12.5	24.4	0.2
Apiary 4 (Chickasaw, NAG)	H13	–	–	3.2	–	–	–	87.0	11.4
	H14	–	–	–	–	–	–	–	–
	H15	0.3	–	–	–	–	16.6	91.6	7.5
	H16	–	–	–	–	7.3	17.6	98.4	9.1
Avg		0	–	0.8	–	1.8	8.5	69.2	7

(–) means pollen not found.

Various noncrop pollens of shrubs and herbaceous plants were identified. The most common noncrop pollen sources are reported in Table 2 and include: 1) Trifolium sp. (clover), 2) Sambucus sp. (elderberry), 3) Rhus typhina (sumac), and 4) Plantago sp. (plantain). Sumac appeared to be the most abundant source of pollen in the NAG area, representing 69.2% of the pollen identified (Table 2). Assorted noncrop pollens were identified at low percentages in Jackson, assumed to originate from urban landscapes.

Pollen Chemical Residues

Sixteen trapped pollen samples (four per apiary) were chemically analyzed (Table 1). Only positive results are reported in Table 3, but full analytical reports are available in the DOI of the supporting information. Low concentrations of fungicides, herbicides, and insecticides were identified at levels well below the oral LD50 values for honey bees (Table 3). Only one pesticide, a fungicide (trifloxystrobin), was detected at 190 PPB in pollen collected from the NAG area (Chickasaw). Imidacloprid was the only neonicotinoid insecticide detected, and it was found (3 PPB) in pollen collected by bees in the most intense AG area (Yum-Yum; Table 3).

Table 3

Results of the pesticide residue detection performed by LC-MS for trapped pollen, of each location. LD50 is based on the data provided by Sanchez-Bayo and Goka (2014) and the toxicity databases ECOTOX of the US Environmental Protection Agency

Sample	Pesticide	Apiary 1	Apiary 2	Apiary 3	Apiary 4	LOD	LD50 Oral (ng/bee)
		Jackson	Milan	Yum-Yum	Chickasaw	PPB
Trapped pollen	Imidacloprida	–	–	3	–	1	13
	Azoxystrobinb	12	88	–	–	2	25,000
	Pendimethalinc	–	57	–	–	6	665,000
	Oxamyld	–	20	–	–	5	380
	Carbendazimb	–	–	5	–	5	50,000
	Atrazinec	12	21	7	–	6	1,000
	Trifloxystrobinb	–	–	–	190	1	200,000
	Chlorpyrifosd	72	–	–	–	1	130

Insecticide.

Funigicide.

Herbicide.

Acaricide.

(–) means no pesticide detected.

Generalized Mixed Effects

Eight different GLM models were built to test the colony weight and brood production as a function of the location, temperature, brood and weight depending on the variable response used in each model (Table 4). Regardless the apiary location, both variables (weight and brood) are correlated and significantly affect each other. Model 7 and 8 estimates indicate respectively that an increase of one brood cell leads to a (0.0005 kg) gain in colony weight, while a gain of (1 kg) in colony weight increases the brood production by 224 cells (P < 0.001; Table 4).

Table 4

Output results of the generalized linear mixed-effects analyses conducted on the dataset. (X ∼ Y1 + Y2…) Respond and explanatory variables as well as the fixed effects used in every model are reported. In all GLMMs, one random effect was considered, which is the apiary being in AG or NAG areas

Model No./ Response ∼ Explanatory variable	Fixed effects	Estimate value	SE	df	T-value	P-value
1) Weight ∼ location	Yum-Yum	6.11	1.7	285	3.5	< 0.001***
	Milan	3.71	1.3	285	0.4	0.008 **
	Jackson	0.68	1.4	285	5.3	0.62
2) Weight ∼ location+brood	Brood	0.0005	0.0001	87	3.8	< 0.001***
	Yum-Yum	3.69	1.78	87	2	0.04 *
	Milan	0.52	1.73	87	0.3	0.76
	Jackson	−0.8	1.72	87	−0.4	0.67
3) Weight ∼ location+brood+temp	Yum-Yum	3.7	1.8	70	2.6	0.01*
	Brood	0.0004	0.0001	70	2.6	0.01*
	Temperature	0.44	0.34	70	1.2	0.2
	Milan	0.72	1.9	70	0.3	0.7
	Jackson	−0.98	1.9	70	−0.5	0.6
4) Brood ∼ location	Jackson	3706	1039	104	3.5	< 0.001***
	Milan	3302	1059	104	3.1	0.002 **
	Yum-Yum	2244	1039	104	2.1	0.03 *
5) Brood ∼ location+weight	Weight	249	64	87	3.8	< 0.001***
	Jackson	3298	1070	87	3	0.002 **
	Milan	2610	1094	87	2.3	0.01 *
	Yum-Yum	413	1187	87	0.3	0.72
6) Brood ∼ location+weight+temp	Jackson	3371	1143	70	2.9	0.004**
	Temperature	546	207	70	2.6	0.01*
	Weight	183	70	70	2.6	0.01*
	Milan	1744	1183	70	1.4	0.1
	Yum-Yum	538	1176	70	0.4	0.6
7) Weight ∼ brood	Brood	0.00053	0.00014	90	3.6	< 0.001 ***
8) Brood ∼ weight	Weight	224	65	90	3.4	< 0.001 ***

The location (landscape) significantly affected colony weight; the model 1 shows that high AG (Yum-Yum) location had the highest impact on colony weight (Estim. = 6.11, P < 0.001), followed by Milan (Estim. = 3.71, P = 0.008), while the Jackson location had neutral effects on colony weight (P = 0.62), Table 4. When brood and temperature components were added to the model 1 (Model 2 and 3), location effects in Jackson still had no influence on the weight variable (P = 0.6). The brood production parameter, however, had higher significant impact on the weight development than location and temperature (Model 2 and 3: brood = P < 0.001 and P = 0.01 respectively), Table 4.

Brood production, as a response variable, was significantly affected by location of colonies. For instance, Jackson location (AG area with urban activity) had the highest levels on brood production (Estim. = 3706, P < 0.001) followed by Milan (P = 0.002) and Yum-yum (P = 0.03), Table 4. Finally, when all variables are included to evaluate the brood production, the only location that positively impacts the brood production is Jackson (Model 6, P = 0.004), with relatively less importance to the temperature and weight in enhancing brood production (P = 0.01, Table 4).

Discussion

Colony biomass, an important but not a perfect proxy to assess colony performance, showed interesting variability within and among the different locations (Fig. 2 and 3). Honey bee colonies located in relatively high AG areas (Milan and Yum-Yum) increased in weight much more than those in the NAG area (Chickasaw, Fig. 3A). The most likely explanation for these differences is better nutrition sources and nectar yields available in the AG areas. This assumption is further supported by the fact that colonies located in Jackson (highest urban habitat and low AG area) did not differ significantly from Chickasaw’s colonies (P = 0.94; Fig. 3A). This result is also confirmed by the GLM analysis, in which Jackson showed no effect (P = 0.62) vis-à-vis the weight gain contrary to Yum-Yum (P < 0.001) and Milan (P = 0.008; Model 1, Table 4). Contrary to our finding, a recent study documented significantly higher colony weight gains in landscapes composed of more than 50% urban areas than colonies situated in areas comprising 50% and more agricultural activity in Denmark (Lecocq et al. 2015). It is, however, difficult to fairly compare both studies, as Lecocq et al. (2015) did not clearly document the nature of the flowering crops as well as differences in the percentages of the agricultural areas in both studies.

The weight correlation matrix of Fig. 3A differentiated three colony clusters: Colonies H(2, 14, 15, 16), H(13, 10, 12, 11), and H(3, 4, 1, 9, 8, 5, 6, 7), colonies of each cluster shared similar weight patterns or behavior. With the exception of H15, all the hives in the first cluster (i.e., weak hives) died during the experiment. The second cluster characterized colonies that exhibited the best long-term weight stability (not necessarily the highest in weight value). The third cluster represented colonies that expressed queenlessness and the highest weight fluctuations (Fig. 2 and 4, Heat map; Fig. 3A). From a location or landscape viewpoint, those clusters correspond to Chickasaw, Yum-Yum and (Milan + Jackson) respectively, with significant differences between Chickasaw and Yum-Yum/Milan colonies (P < 0.001; Fig. 3A). This clustering again confirms the influence of the various landscapes on honey bee performance and that colonies located in the highest AG area had the best long-term weight stability.

Brood production was correlated with total colony biomass. ANOVA tests confirmed that the experimental hives were not rearing brood at similar rates in each location (Fig. 3B). Two AG area apiaries (Jackson and Milan) produced more brood than the NAG area apiary (F = 5.42, n = 128, P = 0.001; Fig. 3B). ANOVA and GLM analyses were in agreement that the Jackson location had the most influence (P < 0.001) on brood production followed by Milan (P = 0.002) then Yum-Yum (P = 0.03; Model 4, Table 4).

It was apparent that bees located in AG areas had access to higher and more sustainable sources of nutrition than those of the NAG area, and starvation losses were only observed in the NAG area. However, mortality of foraging bees resulting from foliar pesticide applications were documented at the Milan and Jackson locations on several occasions (data not shown). These losses did not have measurable impacts at the colony level such as brood production or hive weight. The particular, negative correlations of H2 and H9 with the rest of the colonies are mostly because of divergences in brood production with other colonies as both colonies experienced queenless periods (Fig. 4 and 6B).

Pollen identification clearly indicated that bees intensively foraged or encountered all available crops except cotton (Table 2). This is not surprising, as bees rarely collect cotton pollen (Vaissière and Bradley 1993), but it is unclear how much bees utilized cotton as a source of nectar. The presence of a small amount of sorghum pollen (0.8%) in the NAG area confirmed that bees must have collected pollen from sorghum fields located at 3 km east of Chickasaw State Park (Table 2). Pesticide residues in pollen did not occur at levels expected to cause meaningful mortality to honey bees, and as might be expected there were fewer contaminants found in pollen from the NAG area (Table 3). Similar studies report very weak or insignificant pesticide residues in trapped pollen as well. Alburaki et al. in 2015 documented less than 1 PPB of carbaryl (an insecticide) in their total pollen samples in a Canadian study testing the impact of the pesticides on honey bee health. Similarly, in western Tennessee, only one of 22 pollen samples collected from returning foragers were tested positive for neonicotinoids at trace level (< 1 PPB; Stewart et al. 2014). In a recent comparative assessment study of apiaries in urban, rural, and agricultural areas, the highest level of neonicotinoids found in beebread was 3.9 PPB (imidacloprid; Lawrence et al. 2016). Another study indicated that honey bees foraging in large canola (Brassica napus) fields in southern Ontario (Canada) collected pollen containing 0.5–2 PPB clothianidin (Cutler et al. 2014).

Exposure to pesticides in colonies of the AG areas did not result in measurable impacts on colony productivity. There may be other effects not measured in this study, such as recent evidence that low concentrations of neonicotinoids (4.5 PPB thiamethoxam and 1.5 PPB clothianidin) can reduce honey bee drone reproductive quality (Straub et al. 2016). Indirect impacts of the agricultural pesticides on honey bee performance and health have also been reported (Di Prisco et al. 2013, Alburaki et al. 2015, Di Prisco et al. 2016). Although the NAG environment had fewer pesticide contaminants, bees were challenged to find sustainable food resources, which explains why Chickasaw’s colonies lagged behind other locations in average weight. Two colonies in the NAG area (H13, H16) simultaneously collapsed because of starvation at end of March 2016 (Fig. 6B), which were the largest and most populated colonies during the previous summer. Beekeepers would immediately provide a supplemental sugar source for their starving hives, but this was not done in our case to maintain treatment neutrality. Although not addressed in this research, starvation events may be attributable to poor adaptation of the subspecies to local flora and environment (Louveaux 1973).

The optimal brood nest temperature is known to be close to 35 °C (Seeley and Heinrich 1981, Southwick 1985). Colonies maintained in AG areas had a significantly higher colony temperature (31.33 °C) compared to colonies in the NAG area (29.69 °C). This result is similar to what was obtained for colony weight and brood between AG and NAG hives (Fig. 3A, B). Although the genetic diversity of honey bee colonies are involved in their capacity to accurately regulate temperature and maintaining homeostasis (Jones et al. 2004), the observed variance is more likely related to the difference in population size and activity as our experimental queens were derived from the same genetic source. GLM analyses provided additional interesting results in that, when all variables are included, colony temperature had a higher impact on brood production (Model 6; P < 0.01) than on colony weight (Model 3; P = 0.2), Table 4. Moreover, the same models (3 and 6) indicated significant reciprocal effects between both weight and brood variables (P = 0.01), which is more pronouced in model 7 and 8 when location is disregarded (P < 0.001), Table 4. The Jackson location deserves particular attention regarding brood production. Three GLM models (4, 5 and 6) indicated positive effects of this location on brood productivity (Table 4). One possible explanation for this interesting finding could be related to the elevated urban activity in this location, which provided bees with higher diverse pollen, which was described to enhance colony development (Antunez et al. 2015).

Across locations and the duration of this study, 5 of the 16 colonies were lost (Table 1). Collectively, a 30% colony loss is consistent with loss rates commonly reported by beekeepers (Seitz et al. 2015). The cause of death among colonies varied depending on location, and it is important to interpret each case independently. For instance, the loss of NAG colonies (Chickasaw) was almost certainly not the result of exposure to pesticides, and starvation symptoms were seen and recorded on those colonies. However, other symptoms such as queenlessness, weakened hives, and mosaic brood were the main reasons for colony loss in AG areas.

In conclusion, honey bee colonies foraging in moderate and high AG areas were clearly able to grow faster and to a larger size as a result of better access to sustainable nutrition sources than bees foraging in NAG area and a low AG area with urban activity. However, only the low AG environment with urban activity showed positive effects on brood production when all variables are accounted. Although negative effects of pesticide on colony health were not detected, sublethal doses of insecticides and fungicides were identified in trapped pollen. Better nutrition sources and nectar yields in AG areas helped to develop greater population size, which in turn enabled better colony thermoregulation. NAG areas may provide a less-toxic environment for honey bees but might not provide sustainable foraging resources, leading to colony starvation. Thus, there appears to be a trade-off between increased food resources and the potential for exposure to pesticides in agricultural systems. Careful selection of pesticides and conscientious application of bee-toxic pesticides should greatly reduce the risk of honey bee exposure. From a pollinator enhancement perspective, some non-crop flowering plants were identified as bee-attractive plants and can be enhanced particularly in noncrop areas.