Landscape factors influencing honey bee colony behavior in Southern California commercial apiaries

Colony brood levels, frames of bees (adult bee mass) and internal hive temperature were monitored for 60 colonies for each of two years as they were moved from agricultural, tree crop and mountain landscapes in southern California to blueberry and almond pollination sites. Hive weight was also continuously monitored for 20 of those hives for 6 weeks for both years, during commercial pollination. Pesticide residues in wax, honey and beebread samples were analyzed by composite apiary samples. While colonies in mountain sites had more adult bees and brood than those in agricultural sites in August, by October brood levels were higher in colonies from agricultural sites. Though hives from different original landscapes differed in size in October, hive assessments revealed no differences between the groups after co-wintering when graded for commercial almond pollination. Beebread from hives in agricultural sites had greater agrochemical diversity and in general higher pesticide hazard quotients than those from mountain sites, but those hives also had higher and more constant temperatures from September until January than hives from mountain sites. Hives placed in commercial almond pollination gained on average 287 g per d, compared to an average loss of 68 g per d for colonies in commercial blueberry pollination, although weight data indicated greater foraging effort by colonies in blueberries, possibly due to the proximity and abundance of almond pollen during bloom. Temperature monitoring was effective at distinguishing hive groups and had the best overall value in terms of equipment, installation, colony disturbance and information yield.

www.nature.com/scientificreports www.nature.com/scientificreports/ sites had no significant commercial agriculture with foraging distance; main forage plants, observed by the team were California buckwheat (Eriogonum fasciculatum), mountain sage (Salvia regla), and broomweed (Gutierrezia sarothrae) (D. Winter pers. comm.). MT apiaries the second year were slightly higher (1100 m) than those the previous year (800-900 m) and forage included redshank (Adenostoma sparsifolium) as well as buckwheat, mountain sage and broomweed (D. Winter pers. comm.).  Table 1. Average percentage (±standard error) for a circle with a radius of approximately 1.7 km (=approximately 910 ha) of land around the apiaries in this study (n = 4 for AG, CA and MT; n = 1 for blueberries and almonds), according to the Cropscape web site. "Blue" = blueberry farm, "Alm" = almond orchard, "Dev. " = developed, "intens. " = intensity. Values without standard error indicate that only one sample had the reported use category.
www.nature.com/scientificreports www.nature.com/scientificreports/ In summary: 1. AG sites were about 78 ± 11% field crops, <1% tree crops and 22 ± 11% other kinds of land use; 2. CA sites were about 6 ± 0% field crops, 23 ± 2% tree crops and 77 ± 2% other kinds of land use; and 3. MT sites were <1% field or tree crops and almost entirely other kinds of land use. Some caveats are in order with respect to the analysis using Cropscape data, at least in the case of the "blueberry farm" near Valley Center. Cropscape did not indicate any blueberries at that site, although blueberry production is a layer on the Cropscape map, but numerous visits by the research team to that area indicated several ha of blueberries.
Response variables with respect to bee hives in this study included frames of bees (FOB), brood surface area (cm 2 ), Varroa density (mites per 100 bees), hive temperature, divided into average temperature and temperature variability. In addition, continuous weight data were gathered on 40 hives, and that data used to estimate average daily weight change, nightly hive weight change, the start (dawn) and stop (dusk) of the daily colony activity cycle, and foraging effort. Apiary-level data included land use estimates and agrochemical diversity and concentrations in hive products (wax, honey and bee bread). All raw data on land use categories, hive assessment, temperature and weight are available in the Supplementary File. Hive assessment. Brood surface area changed across the three time points (August, October and March) ( Table 2 and Supplementary Table S1). The three original landscapes, Imperial Valley agriculture sites (AG), sites with comparatively large acreage of citrus and avocado orchards (CA) and non-agriculture higher altitude sites (MT) were significantly different in August, at the start of both experiments. The AG colonies had significantly less brood surface area than either the CA or MT colonies. In October brood surface areas were again significantly different, but the AG colonies had significantly more brood than the MT colonies. Both August and October brood production were also significantly different with respect to Year. On the last sample date in March, after almond pollination, no significant differences were observed among the original landscapes.
Comparisons of frames of bees (FOB) showed differences among original landscape groups, but only for the August inspection (Table 1 and Supplementary Table S2). As with brood surface area, colonies in the MT group had significantly more bees than did colonies in the AG group. FOB data could be gathered without serious hive disturbance, so data were gathered in January as well, and the October to January data were included in the same analysis. Neither original landscape nor year were significant in the October-January analysis. No fixed factors were significant after almond pollination.
Hive temperature. Within-day temperature amplitudes from the beginning of the experiment in mid-August until the end of October were significantly affected by original landscape across both years of the study but average daily temperature was not (Figs. 2 and 3, Supplementary Table S3). Post hoc contrasts showed that bee colonies kept in the Imperial Valley during the summer had lower temperature amplitudes during that period than colonies placed in the CA environment (Table 1). Restricting the analysis to November to mid-January, when the challenge to maintain temperatures within the hive is greatest, the AG colonies maintained average temperatures 3.9 to 4.2 °C higher, and temperature variability 0.7-1.4 °C lower, than colonies at either of the other sites over both years of the study. Considering the entire period from mid-August to mid-January and using the "slices" function in SAS, significant treatment effects were observed with respect to average hive www.nature.com/scientificreports www.nature.com/scientificreports/  www.nature.com/scientificreports www.nature.com/scientificreports/ temperature for 74 of the next 80 days after 29 October, but not before that date. Significant treatment effects were observed with respect to temperature amplitude for 112 of 115 days after 16 September.
Neither temperature average nor amplitude were significantly different among original landscape groups when analyses were restricted from mid-January to the end of almond pollination in early March, although temperature averages were different between pollination sites (Supplementary Table S4). However, those temperature differences were inconsistent between the two years with temperature averages higher in almonds in 2017 but lower in 2018 ( Table 3). Colonies that failed during the course of the study generally failed in late fall and early winter, and exhibited distinctive temperature patterns when compared to overall colony averages without those colonies (Fig. 4). Colonies were removed from the study once their death had been confirmed by visual assessment.
Hive weight. Daily hive weight changes were significantly different between colonies sent to almond pollination and those retained in blueberry pollination (Table 4, Fig. 5, Supplementary Tables S5 and S6). Colonies placed in almond pollination gained on average 287 g per d between the two years of the study, and those retained in blueberry pollination lost on average 73 g per d. Distinct patterns in increasing and decreasing rates of hive weight change were compared to rainfall data for each location. Rainfall near Escondido was higher and of the seven peaks in hive weight gain across both years, five were associated with the five highest local rainfall events, suggesting that without rainfall, hive weight loss would have been higher. Of nine peaks on hive weight gain in almonds between the two years, only four were associated with rainfall and in almost all cases only partially. The hive weight change patterns in almonds may have been be related to factors such as the variety or local environment of the plants upon which the bees were foraging.
Dawn (t D ), dusk, the night slope (S N ) and the hive weight change due to initial forager departure (ΔF) were significantly different between pollination environments. Dawn break points were roughly consistent between the two years, with the colonies in blueberries beginning activity 20-33 minutes before those in almonds, and no differences were detected between years. Dusk break points were also different, with colonies in blueberries ceasing foraging 7-48 minutes earlier than those in almonds, but dusk times were also significantly different between years.
Night slopes were significantly different between colonies placed in almonds and those placed in blueberries. Hives placed in almonds lost, on average, about 12.6 g/h during the night and those in blueberries lost about 6.7 g/h. These differences likely reflect the quantity of nectar and pollen in the hives; the higher value for hives in almonds simply reflects the larger surface area of drying nectar and pollen than those in blueberries, which was reflected in the higher average daily hive weight change for hives in almonds discussed above. The initial forager mass, after subtracting an estimated value for the moisture loss during that period, were significantly lower for colonies placed in almonds compared to those placed in blueberries, after controlling for adult bee mass.
Mite density. Mite levels were not significantly different among the original landscapes but they were significantly different between years (Supplementary Table S7). Mites were lower in the 2016-17 trial (0.33 ± 0.19 and 0.34 ± 0.24 mites per 100 bees for the August and October samples, respectively) than in the 2017-18 trial (1.94 ± 0.27 and 1.87 ± 0.32 mites per 100 bees for the August and October samples, respectively). After almond pollination in 2017 no mites were detected whereas in 2018 colonies placed in almonds had 0.21 ± 0.08 mites per 100 bees and those in blueberries had 0.73 ± 0.38 mites per 100 bees.  www.nature.com/scientificreports www.nature.com/scientificreports/ Agrochemical residues. Three matrices were evaluated for residues: honey (or nectar), wax and beebread. A total of 68 compounds were detected, including fungicides, herbicides, insecticides, miticides, two synergist compounds and an insect repellent (Supplementary Table S8). Honey samples were evaluated over four sampling occasions in the 2016-17 experiment (Supplementary Table S9). DMPF (a breakdown product of amitraz, a miticide) was detected in all samples at concentrations ranging from 6 to 155 ppb. Flonicamid, a pyridine compound that affects insect chordotonal organs, was also detected in samples only from hives kept in the Imperial Valley and at concentrations ranging from "trace" (<7 ppb) to 20 ppb. The following year only honey samples from the August sampling date were analyzed (Supplementary Table S10). DMPF was again found in all samples. Flonicamid was found in hives from citrus-avocado orchards in near Valley Center, CA, and one sample from the higher altitude hives but not in the AG hives. Methoxyfenozide, an insect growth regulator, was also detected in one of the CA samples.
Wax samples in August also had a high compound diversity (Table 5) as well as very high concentrations of certain compounds, particularly, but not exclusively, miticides. One wax sample from August 2016 had a DMPF concentration of 39,900 ppb, which was the highest concentration of any compound in any matrix in this study   www.nature.com/scientificreports www.nature.com/scientificreports/ (Supplementary Table S11). In those same samples, fluvalinate ranged from 450 to 3720 ppb and coumaphos from 69-370 ppb. Two insecticides, fenpyroximate (up to 98 ppb) and phenothrin (up to 1900 ppb), were both found in all wax samples from that date. Permethrin was found in several samples across all original landscapes at up to 2430 ppb. Trifluralin, an herbicide, was found in most samples at up to 189 ppb. In the August 2017 samples fenpyroximate (up to 98 ppb) and permethrin (up to 1650 ppb) were found in all samples (Supplementary Table S12). Analysis of agrochemical diversity showed that there was no interaction between year and original landscape (P = 0.096) and in the main effects model that original landscape was not significant (P = 0.096).
August, October and March beebread samples in both years were analyzed. The number of compounds detected in beebread prior to almond pollination was high compared to honey (Table 4). More insecticides were detected in the AG hives than in samples from any other sites, with the exception of one sample from CA in 2017-18, and more fungicides and herbicides were detected in CA and MT hives than those from the AG hives. Methoxyfenozide was detected in the most samples and at the highest concentration (350 ppb) in 2016-17 (Supplementary Table S13). The August 2017-18 sample from the MT 4 site had very high levels (2220 ppb) of cyphenothrin, a pyrethroid insecticide (Supplementary Table S14), and samples from the AG hives on that date had high methoxyfenozide concentrations, up to 1520 ppb, as that original landscape did the previous year. Also in 2017-18 both MT and AG samples had low levels (4-7 ppb) of diethyltoluamide (DEET, an insect repellant). Samples from October 2017 showed high levels of the fungicide chlorothalonil at 398-908 ppb among the AG3, AG4, MT3 and MT4 sites, with none reported from CA3 or CA4 (Supplementary Table S15). The MT3 sample also had 536 ppb of prodiamine, an herbicide not found in any other samples. An ANCOVA analysis of agrochemical diversity for the August samples, with "year" as the covariate, showed that original landscape was significant (P = 0.039). Post hoc contrasts indicated that the AG sites had significantly more kinds of compounds (13.3 ± 1.4) than the MT sites (7 ± 0.7) (no other contrasts were significant).
Samples from colonies sent to almonds had both a high diversity of compounds as well as high concentrations of certain compounds (Supplementary Table S16). Samples from March 2017 had pyriproxyfen (an insect growth regulator) at 480 and 870 ppb, as well as chlorpyrifos at 31-83 ppb and the fungicide cyprodinil at up to 48 ppb. Samples from almonds in March 2018 had the fungicide chlorothalonil at 250 to 6860 ppb and the herbicide pendimethalin at 110 to 132 ppb, as well as chlopyrifos in one sample at 287 ppb (Supplementary Table S17). With respect to agrochemical diversity for the March samples, there was a significant interaction between year and original landscape (P < 0.001); main effects of year and original landscape were also significant (P < 0.001 for both). A higher diversity of agrochemicals was observed in 2018 compared to 2017, and in almonds compared to blueberries (Table 6).

Discussion
In this study commercial honey bee colonies spent several months in apiaries in an "original landscape", after which they were moved to a common site. If the original landscape had a short, limited impact on colony growth and behavior, then colonies would be expected to exhibit roughly the same behaviors among themselves after having been moved to the common landscape, or shortly thereafter, irrespective of the original landscape. After www.nature.com/scientificreports www.nature.com/scientificreports/ several months in the common site, half the colonies were moved to almond pollination and half remained at the common blueberry pollination site. Our goal was to monitor colonies during a typical commercial migration among sites and observe how colonies reacted and adapted to these changes.
The FOB visual assessment detected landscape-level differences among colonies only for the first sampling occasion, in August, across the two years of the study. Colonies from MT sites were, on average, about 10% larger than those from the AG group based on the FOB data. Brood levels were significantly higher among colonies from MT sites than those from the AG sites in August but the opposite was true by October. Average within-hive temperatures were slightly but consistently higher, and temperature amplitudes lower, in the AG colonies than either the CA or the MT colonies until mid-January. However, significant landscape effects on temperature amplitude were only detected starting mid-September (and continuing, with a few exceptions, for each time point until January) and temperature averages starting the end of October until January.
By mid-March no differences were observed between colonies in blueberries and those in almonds pollination with respect to either FOB or brood levels. Average hive temperatures were significantly different between the two pollination sites, but the direction of the difference changed between the two years, suggesting that it may have been due to weather or other environmental factors. The list of potential factors that might give rise to observed differences among the original landscape groups was large and investigative resources limited, so the study focused on agrochemical residues in honey, beebread and wax, and Varroa mite density.  Table 5. Number of agricultural compounds by type reported from wax and beebread for 2 apiaries per original landscape across 2 years, and beebread hazard quotients. In October 2016 samples were pooled between apiaries within original landscape so results for just that level is reported. "1 st " and "2 nd " refer to the two apiaries within a given landscape for a given year. The synergists piperonyl butoxide and MGK264 were included as "insecticides" for counting purposes.
www.nature.com/scientificreports www.nature.com/scientificreports/ Hive weight data were consistent between years, with hives in almonds gaining weight both years and those in blueberries losing weight in both years. Both the start and end of daily activity were significantly earlier for colonies in almonds compared to those in blueberries, but that was likely related to physical characteristics of the apiaries and how the hives were placed. The direction a hive is facing can affect the timing of the daily activity period 25 . Much of the hive weight gain during these periods was probably due to nectar and pollen collection. Hive weight loss at night during a nectar and pollen flow has been attributed to moisture loss of the drying nectar and pollen 22 . Hives in almonds lost on average more than 12 g per h at night while hives in blueberries lost on average less than 7 g per h, which suggests that, assuming the water content of the collected nectar was similar between sites, there was more material to dry in the hives in almonds.
Initial foraging populations were estimated by subtracting the estimated weight loss due to drying from the total hive weight change during the period of initial forager departure in the morning. These forager populations were significantly smaller for colonies placed in almonds compared to those in the blueberry environment, in spite of the higher weight gain for those colonies, which may have been due to a number of factors. One factor is the high density of forage in the almond environment -the colonies in almonds may have been able to maintain a higher foraging success with a smaller forager population compared to colonies in blueberries because of the high density of forage in almonds and/or because almond pollen is highly preferred compared to pollen in the blueberry site. Also, these estimates were made from changes in hive mass over time. It is also possible that forage was much more proximal in almonds, and surfeited foragers started returning as other foragers were still making their initial departure, confounding departures with returns and resulting in an underestimate of the forager population.
Mite levels were low overall and not significantly different among original landscapes, indicating that differences in mite densities were unlikely to explain differences among treatment groups. Varroa mite infestation is strongly correlated with some important viral diseases 26 -we assumed a high incidence of at least Varroa-transmitted viral diseases was less likely when mite levels are as low as they were observed here. However, other diseases could also have been differentially distributed.
Pesticide residues were measured as composite samples on the level of the apiary or landscape. Honey samples were found to have few compounds -mainly DMPF, a by-product of the miticide amitraz, and flonicamid, an insecticide. Comb wax, which was only analyzed for the August samples in both years, typically had several compounds at comparatively high levels (>500 ppb). However, the relationship of residue concentrations in wax and bee health has not been firmly established, at least for many of the compounds detected. Wax, being lipophilic, can trap and store compounds, which may remain for years with little exposure to light, moisture or microbial activity. Because the compounds may be embedded in the wax matrix, the exposure of bees to those compounds may also be limited. Averaging across original landscapes and years, colonies kept in the AG and CA sites had on average 15.5 compounds per wax sample, while those in MT sites had on average 12.8 compounds, which is consistent with land use patterns.
Beebread may be the most revealing matrix for estimating the concentrations and diversity of agrochemicals bees are exposed to at a given point in time 6 . A recent study found that among 6 geographically disparate apiaries in the US, 79 different pesticides and metabolites were observed, with up to 10 distinct modes of action 27 . While wax and honey can last for many months or even years, bees prefer pollen <72 h old 28 , so it would more likely reflect the hive's recent environment than wax or honey. In both August and October, beebread samples from the AG sites had more compounds (13.3 and 11.8, respectively) than either those in CA (11 and 5.8, respectively) or MT landscapes (7 and 8, respectively). These compounds included neonicotinoids, such as imidacloprid, and insect growth regulators (IGRs), such as methoxyfenozide, both of which have been shown to affect thermoregulation in bee colonies 8,25 . Some compounds, such as methoxyfenozide and flupyradifurone, were at high, albeit sublethal levels 29 . The IGRs methoxyfenozide and pyriproxyfen have been shown to reduce forager survival 4 . By October fewer pesticides were being detected overall, but the hives in MT landscapes were clearly being exposed  www.nature.com/scientificreports www.nature.com/scientificreports/ to high levels of fungicides (chlorothalonil), herbicides (chlorthal-dimethyl and prodiamine) and IGRs (methoxyfenozide). Bee colonies were exposed to considerably more different compounds in almond pollination (16.3

, averaging across all landscapes and across both years) than in blueberry pollination (5.3) at that time of year.
Agrochemical diversity is an important aspect of the in-hive pesticide exposome -greater compound diversity increases the possibility of interactions among compounds, which may be significant 6,30,31 . However, compound diversity alone is not likely to be a good measure of risk to bees. The pesticide hazard quotient was developed to take into account compound toxicity and residue concentration 5,6 . Hazard quotients calculated here showed high values (values >1000 are considered "high" 6 ) for colonies in agricultural areas in the first year, but lower quotients for the same locations in the second year. In October quotients were low overall. Two mechanisms for lower agrochemical concentrations in the composite samples are: (1) bee colonies had largely consumed the pollen collected in the original landscape, given their preference for freshly-collected pollen, and replaced it with pollen from the new site; and (2) agrochemicals had decomposed perhaps through microbial activity. The second mechanism is unlikely, as microbial cell densities in bee bread are generally very low and static, even decreasing over time 32 . Some samples from the colonies kept in non-agricultural areas (MT) showed very high quotients due to high levels of particular compounds, suggesting the bees in at least some colonies may be been subjected to lethal levels. For samples collected in March, the composite bee bread sample from colonies used for almond pollination had high levels the first year, particularly compared to samples from colonies that remained in blueberries, but much less so the second year. No unusual bee mortality, such as masses of dead adult bees, was noted in any yards at any time.
Another approach would be to consider where and when compounds of known toxicity to honey bees were detected (Table 7). Emamectin benzoate and pyridaben are both comparatively toxic to honey bees but both those compounds were only detected in samples collected in March 2017 (Emamectin benzoate from colonies in blueberries and pyridaben from colonies in almonds) and thus their distribution has little explanatory value. Permethrin, found only in wax samples, is also somewhat toxic to bees 29,33 but it was reported from at least one sample from all original landscapes in both years, so its occurrence would not explain landscape differences. A fourth example, chlorpyrifos, is comparatively toxic to bees 29,33 and has been found to affect learning behavior at sublethal concentrations 34 . Chlorpyrifos was detected in beebread in August (both years) only from colonies in the AG group, and again in samples taken in March only from colonies used for almond pollination. Colonies from the AG sites had the highest average brood temperatures and lowest temperature variability, and colonies from almond pollination had the highest daily weight gain, both of which are generally considered positively correlated with colony health.
The observed relationship between chlorpyrifos exposure and colony thermoregulation and weight gain suggests several non-exclusive possibilities, three of which are listed here. Some colony effects may be due to hormesis, defined as a change in the shape of the dose-response curve at low, sublethal concentrations of toxic compounds 35 . Compounds may cause negative effects at higher concentrations but not at low concentrations. Another possibility is that the LD 50 may also be a poor measure of the impact of a compound on colony-level behavior. Several researchers have explored the use of hazard 36 or risk quotients 37 that take into account acute toxicity as well as time of exposure, insect life stage most at risk of exposure, application method and other factors (e.g., https://www.epa.gov/pesticide-science-and-assessing-pesticide-risks/models-pesticide-risk-assessment). Also, some kinds of exposure may affect honey bee behavior without affecting survivorship. Imidacloprid, which has a comparatively low LD 50 25,33 , had measurable effects on colony behavior at concentrations that had no apparent effect on adult bee longevity 8 raising the question of whether there are compounds with higher LD 50 values but that also affect colony behavior. Incorporating colony-level behavior studies when evaluating novel pesticides may need to be considered.
A third possible explanation is that the observed distribution of vigorous colonies was not due to agrochemical exposure but rather to other factors such as ambient temperature or food quantity and/or quality. Colonies in the AG sites had greatest access to alfalfa fields; exposure to alfalfa has been associated with lower colony strength in the long term but positive colony growth in the short term 36 and those colonies would also have had access to non-crop plants, particularly along irrigation canals and drainage ditches. Even together on the blueberry farm, colonies may still have had variable diets. Bee colonies in commercial blueberry fields in Canada have been found to have a low proportion of crop pollen compared to non-crop pollen and that beebread from those colonies has relatively low nutritional value compared to, for example, colonies placed in apple orchards 30 . If the bees did not prefer blueberry pollen they may have sought forage elsewhere. Finally, colonies in almond pollination gained weight while those in blueberry fields lost weight. Bee colonies in particular agricultural landscapes may thus thrive in spite of increased agrochemical exposure.
This work concerns a longitudinal study embedded in a commercial beekeeping operation, and was intended to monitor colony growth and activity in ways that would reveal differences among groups that strictly visual assessments may not detect. Such differences were observed and in the case of thermoregulation, those differences were largely consistent between the two years of the study. Colonies clearly adapt many aspects of their behavior and phenology to their current situation, but the effects of previous exposure to other environments can have be longer term and difficult to observe using strictly visual assessments. Owing to resource limitations (mainly sample processing costs), insufficient data was available to conclusively determine the cause of those differences among landscapes. Per hive analyses of agrochemicals and beebread nutritional quality, e.g., amino acid and essential fatty acid concentrations, and bee disease loads (especially for diseases not associated with Varroa) would have been very informative.
Which data in this study were the most valuable, and which data were the best value, taking into account the cost in terms of time and money resources? With respect to measures of colony health, estimating frames of bees was the least expensive, taking usually less than a minute and little cost other than transportation, but was also the least precise and revealed few differences among treatments. Any differences would probably have to be large www.nature.com/scientificreports www.nature.com/scientificreports/ before they could be detected. Measurement of brood surface area using frame photographs, taken during the visual assessment, had a comparatively low monetary cost (inexpensive cameras can suffice, and the image analysis software was free) and photographing frames can be fairly rapid, but analysis of the resulting photographs was time consuming. Visiting each of 60 hives once could generate 500-800 photographs (depending on the season) requiring on average a few minutes per photograph, and FOB and brood surface area data were collected with 2-4 months between data points. More frequent visits would increase transportation costs and frequency of colony disturbance. Continuous hive weight data showed significant treatment differences, in colony behavior, during a period when no other response variables did, and weight data is easy to collect without disturbing the colonies, easy to analyze and, for the most part, to interpret. However, the precision hive scales used here were expensive and cumbersome to install; less precise scales may be more convenient but may be less effective at detecting differences. Temperature data revealed many differences over a long period, the sensors are fairly inexpensive, and collection and analysis of the data are not time consuming. Further work is needed to improve understanding and exploitation of temperature data; changes in thermoregulation have been shown to be due to many factors. Clearly just one kind of data is unlikely to be sufficient in monitoring colony health.  www.nature.com/scientificreports www.nature.com/scientificreports/

Materials and Methods
Overall experimental design. In August of each of the two years of this study a total of 60 commercial hives were selected in two apiaries, each containing 60-80 hives, in each of three different landscapes: field crop and forage agriculture (AG) in the Imperial Valley (Holtville, CA), citrus and avocado groves (CA) in mid-altitude areas (Escondido and Valley Center, CA) and natural forage at higher elevations (MT) (near Santa Ysabel and Boulevard, CA) ( Table 8). Apiary sites were assessed using the National Agricultural Statistical Service (NASS) Cropscape web site (https://nassgeodata.gmu.edu/ CropScape) to obtain acreage estimates for all land use categories within an approximately 1.7 km radius of the apiary (see Supplementary File). The hives, each of which had two "deep" Langstroth boxes (approximately 40.6 × 50.5 × 24.4 cm), had been either assembled the previous May and provided with a new queen (using open-mated queens from a commercial source) or kept as full colonies from the previous year. All colonies in both years of the study were in their respective sites by the end of May.
Ten hives were selected in each of two apiaries at least 2 km apart within each landscape. Hives were at least 3 m apart within the original apiary and did not neighbor any other selected hives. Frames of bees (FOB) were visually estimated in both top and bottom boxes by the number of between-frame spaces and half spaces covered by bees when observed from above. The same member of the team (EB) estimated all FOB values. The same worker estimated FOB levels in another study in which adult bee masses were measured independently using another method (see 25 ), and the linear regression between that worker's FOB estimates and adult bee masses are provided (see Supplementary File). After visual estimation, both sides of all frames with brood were photographed using a 16.3 megapixel digital camera (Canon Rebel SL1, Canon USA, Inc., Melville, NY). The area of sealed brood per frame was estimated later from the frame photographs using either ImageJ software (version 1.47. W. Rasband, National Institutes of Health, USA), or CombCount; 38 this method has been described in other publications (see 8,13,16,25 ). Any colony that had capped brood was included in the study.
Samples of honey (approx. 2 ml per hive), beebread (from approx. 15 cells per hive) and comb wax (approx. 1 g) were collected. Adult bees (50-100) were collected from the outer brood nest (usually the first frame encountered with brood when working from one side of the box to the other) and stored on ice. Based on where the sample was collected, it was assumed to contain a high proportion of nurse bees with some foragers. A temperature sensor (iButton Thermochron, precision ±0.06 °C) enclosed in plastic tissue embedding cassettes (Thermo Fisher Scientific, Waltham, MA) was stapled to the center of the top bar on the 5th frame in the bottom box of each hive and set to record every 30 min. In October all hives were moved to a blueberry farm (Valley Center, CA) for commercial pollination. During September and early October all hives were given one or more supplemental soybean flour-based protein patties and a 50 g patty consisting of palm oil and 0.09% (w/w) amitraz to control Varroa mites.
In October all hives were assessed again in the same manner, and the temperature sensor was removed, replaced with a fresh sensor, and then downloaded. In December and January 2017 all hives were given protein patty and another patty containing amitraz. In January 2017 all hives were visited to 1) estimate the frames of bees; and 2) change the temperature sensor. Because photographing brood frames is an invasive procedure that requires removing the frames and shaking them free of bees, it was not conducted in January. Hives were randomly assigned to either remain in the blueberry farm or be sent to an almond orchard for commercial pollination. Because commercial pollination contracts stipulate the required size of the bee colonies (frames of bees as estimated by either the commercial beekeeper or another party), not all hive designated for almond pollination were sent to almond orchards. However, all hives designated to remain in the blueberry farm remained there through March. The apiary in the almond orchard contained only hives used in this study (about 30 hives).  www.nature.com/scientificreports www.nature.com/scientificreports/ During the almond pollination period, ten hives were selected from the group of hives sent to the almond orchard and ten selected from the group of hives that remained in the blueberries. Hives were placed on stainless steel electronic scales (Tekfa model B-2418 and Avery Weigh-Tronix model BSAO1824-200) (max. capacity: 100 kg, precision: ±20 g; operating temperature: −30 °C to 70 °C) and linked to 12-bit dataloggers (Hobo U-12 External Channel datalogger, Onset Computer Corporation, Bourne, MA, USA) with weight recorded every 15 minutes. The scales were powered by deep-cycle batteries connected to solar panels. All surviving hives in both groups were evaluated for a final time after the end of the almond pollination period.
Wax, beebread and honey samples from the August evaluation were pooled by apiary (i.e. 6 groups per year) and submitted to the Laboratory Approval and Testing Division, Agricultural Marketing Service, USDA (LATD), Gastonia, NC, for full panel testing (the number of compounds detected during the full panel testing varied from 200-214 compounds, including insecticides, fungicides, miticides, herbicides, synergists and breakdown products; please see Supplementary Tables for a complete list). Adult bee samples were stored at −20 °C and subjected to an alcohol wash to remove mites. Mites and bees were then counted to estimate the number of mites per 100 bees for each sampling date. Internal hive temperature data were divided into daily average and within-day detrended data. Detrended data were calculated as the difference between the 25 hour running average and the raw data. Sine curves were fit, via least sums of squares (see 16 ), to 3-day subsamples of detrended data taken sequentially by day, and curve amplitudes, representing estimates of daily hive temperature variability, were used as response variables. Repeated measures MANOVA was used to evaluate the effects of treatment, day, and their interaction, with the sealed brood surface area at the start of the experiment in August (both years) as a covariate on both the average daily hive temperature and the amplitudes of the fit sine curves. Thermoregulation is closely linked to brood production 16,17 so it was used as a covariate to control for pre-existing differences among colonies. Temperature amplitude datasets were reduced to one value every 3 d for repeated measures analysis to ensure no overlap between subsamples. Temperature data were evaluated with respect to three time periods: mid-August to October, November to mid-January, and mid-January to beginning of March. Colonies were removed from the study once their death had been confirmed by visual assessment.
Continuous hive weight data were likewise considered with respect to daily weight change and within-day changes. Because of the low number of hives on scales per pollination site per year (10), original landscape groups were no longer considered. Daily weight change was calculated by averaging, for each 24 hour period, all the weight data for a given hive and subtracting that value from the value for the following day. These between-day weight changes were analyzed using repeated-measures MANOVA with pollination sites (almond or blueberry; hereafter "pollination"), year, and 2-way interactions as fixed effects and pre-pollination weight as a covariate.
Weight data were detrended for each day by subtracting the hive weight estimate at midnight (or closest time thereafter) from each subsequent weight value over the next 24 h (see 22 ). The resulting within-day weight datasets were modeled using the "segmented" function in R (version 3.6.1. The R Foundation for Statistical Computing) which fits a segmented line derived from a linear or generalized linear model to a dependent variable using a bootstrapping procedure 39 . Piecewise regressions with 4 breakpoints were fit to the data, which yielded estimates for 10 parameters: 4 break point values, 5 slope values and the adjusted r 2 . Because the data were detrended by subtracting the raw data value at midnight, daily datasets were mathematically independent from each other. The following values were obtained: 1. Dawn break point (start of daily hive activity, usually the 1st break point); 2. Dusk break point (end of daily hive activity, usually the 4 th break point) 3. Slope of the first segment after dawn, usually the 2nd segment; 4. Night slope: the slope of the 1st segment, representing hive weight change during the night.
While daily patterns in hive weight change are generally consistent, they can also be variable. To reduce what were considered known errors, if the 1st break point occurred before 4 AM, it was assumed to be an error, the 2nd break point was used instead as the dawn break point (no restrictions were placed on that new estimate) and the 3rd, rather than 2nd, segment was assumed to be the first segment after dawn (value 3 in the list above). Likewise, if the 4th break point occurred after 8 PM it was assumed to be an error and the 3rd break point was taken as the dusk break point (again with no restrictions placed on that new estimate). The slope of the first segment after dawn was multiplied by the length of time between the break points on either end to estimate hive weight loss during that period. This hive weight loss was attributed to two factors: (1) moisture loss (nectar drying, respiration, etc.); and (2) forager departure. Weight loss due to forager departure was then calculated by (1) multiplying the slope of the first segment after dawn by the length of the projection of the segment on the X axis (i.e. time) to generate an estimate of total hive weight loss during that period, and then subtracting from that value the estimated weight loss due to drying: www.nature.com/scientificreports www.nature.com/scientificreports/ where ΔF is the hive weight change due to forager departure, S F is the slope of the piecewise regression segment associated with forager departure, t D is the time of the dawn break point, t D+1 is the time of the following break point, and S N is the slope of the first segment of the regression, representing hive moisture loss due to drying at night (this was used as a proxy to estimate the drying rate during forager departure). Dawn and dusk break points, S N and ΔF were used as response variables in repeated-measures ANOVAs with treatment and year as fixed factors, hive as a random factor and January FOB as a covariate to control for colony size. Days with rainfall >3 mm were removed from the analysis because of the effect of rainfall on the shape of within-day hive weight changes.
Pesticide diversity data (i.e. the total number of compounds detected) in wax and beebread samples were analyzed within sampling occasion using ANCOVA, with year as the covariate and landscape as the main factor. First an interaction model was tested; if the interaction was not significant, slopes were assumed equal and a main effects model was tested. Because beebread samples in October 2016 were pooled across landscape, only August and March beebread data were statistically analyzed. Pesticide hazard quotients were calculated using residue data and information on honey bee LD 50 levels obtained from either the University of Hertfordshire Pesticide Properties Database (https://sitem.herts.ac.uk/aeru/ppdb/en/atoz.htm), or from recent publications, i.e. 6,33,36 . Because "contact" LD 50 data were more readily available than "oral" LD 50 , those data were used throughout. Hazard quotients were calculated by dividing the residue concentration in ppb by the LD 50 value in µg per bee, and then summing the values for all detected compounds for a given sample 5,6 . In cases where the LD 50 was presented as the lower end of a range of values, such as ">100 µg per bee" the inequality sign was simply removed, so the values shown here should be considered lower estimates. Where "trace" amounts of an agrochemical was reported, a token value of 1 ppb was used.