Industrial bees: The impact of apicultural intensification on local disease prevalence

Abstract It is generally thought that the intensification of farming will result in higher disease prevalences, although there is little specific modelling testing this idea. Focussing on honeybees, we build multi‐colony models to inform how “apicultural intensification” is predicted to impact honeybee pathogen epidemiology at the apiary scale. We used both agent‐based and analytical models to show that three linked aspects of apicultural intensification (increased population sizes, changes in population network structure and increased between‐colony transmission) are unlikely to greatly increase disease prevalence in apiaries. Principally this is because even low‐intensity apiculture exhibits high disease prevalence. The greatest impacts of apicultural intensification are found for diseases with relatively low R0 (basic reproduction number), however, such diseases cause little overall disease prevalence and, therefore, the impacts of intensification are minor. Furthermore, the smallest impacts of intensification are for diseases with high R0 values, which we argue are typical of important honeybee diseases. Policy Implications: Our findings contradict the idea that apicultural intensification by crowding honeybee colonies in large, dense apiaries leads to notably higher disease prevalences for established honeybee pathogens. More broadly, our work demonstrates the need for informative models of all agricultural systems and management practices in order to understand the implications of management changes on diseases.

Honeybees are typically managed in apiaries, which are associated colonies placed together for beekeeping convenience at a single site. Pathogen dynamics at the apiary level are determined both by pathogen transmission within and between colonies.
Intensification of apiculture changes apiary ecology in a number of ways, all potentially relevant to disease . In particular, increasing the number of colonies and changing the arrangement of those colonies influences epidemiology through changes in both the size and network structure of the population.
The intensification of agricultural systems generally means larger, denser population sizes and greater pathogen transmissibility at local (within a population, such as a farm) and landscape (between populations, such as neighbouring farms) scales. To understand these effects in honeybees we build multi-colony models to examine how apicultural intensification is predicted to impact honeybee pathogen epidemiology. We examine the epidemiological consequences of increasing the number of colonies within an apiary, changing colony configurations, and increasing between-colony pathogen transmission.

| MATERIAL S AND ME THODS
We combine mathematical models and agent-based model (ABM) simulations to make predictions on how intensification affects disease risk, spread and endemic prevalence within an apiary. The key to our approach is that we capture pathogen transmission both within and between colonies.
We generalize colony arrangements to three unique configurations drawn from experience, classic apicultural literature (Jay, 1966) and current experimental work (Dynes, Berry, Delaplane, Brosi, & de Roode, 2019): array, circular and lattice ( Figure 1). We restrict between-colony pathogen transmission to nearest neighbours (see discussion), those in closest proximity to each other (connected by an arrow in Figure 2). Between-colony transmission is always assumed to be at a lower rate than within colony transmission. The mathematical model allows us to obtain tractable analytical results while the ABM simulations allow us to model disease at the level of the individual bee and consider stochastic effects.
We first derive a compartmental SI (Susceptible, Infected) model for pathogen transmission within an apiary. The model treats each colony as an individual population and allows for within-colony as well as between-colony transmission (for nearest neighbours).
Within a colony, honeybees are either susceptible to infection or infected (and infectious). We denote the number of susceptible honeybees in colony i at time t as S i (t). Likewise, we denote the number of honeybees in colony i infected with the pathogen at time t as I i (t).
Susceptible honeybees in colony i become infected at rate β ij following contact with an infected bee that resides in colony j. We assume that honeybees do not recover from infection. Honeybees are born at rate ϕ, have a natural mortality rate of m and an additional mortality rate of v if infected. The following 2n differential Equations, (1), model disease transmission within and between n colonies in an apiary.
The matrix β = [β ij ] will depend on the colony arrangement (see Figure 1; and S.I. Section 1). The transmission rate between a susceptible and infected honeybee within the colony is a, and transmission between neighbouring colonies is b. For example, for a nine-colony apiary, the transmission matrices for an array, circular and lattice configured apiary (respectively) are as follows: The corresponding network structures for the above transmission matrices can be seen in Figure S1. We assume that honeybees are much more likely to become infected by a honeybee that resides within its home colony than by a honeybee from a neighbouring colony (i.e. a » b). Note that for each apiary configuration to be possible and unique, the number of colonies (n) must be a perfect square, n = L 2 where L ≥ 3 (see Figure 1). Therefore, the minimum number of colonies per apiary is nine, which has been observed to be the mean size of a hobbyist or small beekeeping operation (Mõtus et al., 2016;Pocol, Marghitas, & Popa, 2012).
We can understand the dynamics presented by our models by focussing on the basic reproduction number, R 0 . R 0 is a fundamental concept in infectious disease ecology, defined as the average number of secondary infections caused by one infectious individual in an otherwise entirely susceptible population (Anderson & May, 1992).
We derive R 0 expressions, using model (1), for each of the apiary configurations. R 0 derivations using model (1) allow us to characterize the relationship between R 0 and pathogen prevalence, defined as the proportion of honeybees within an apiary that are infected at the endemic equilibrium. The R 0 expressions for apiaries with n > 1 For the ABM we estimate R 0 values for particular parameter combinations by treating simulation outputs as ideal empirical data (Keeling & Rohani, 2008) and track the number of infections following the index case. The term "base R 0 " is used throughout the remainder of this paper and refers to a value of R 0 for a specific pathogen phenotype in a least intensified apiary, an array with nine colonies (see Figure 2). We determine how intensification affects R 0 by separating R 0 into a "base R 0 " and an "additional R 0 ".
The term "additional R 0 " refers to the observed difference in R 0 for a given pathogen phenotype when comparing a "lower intensity" apiary to a "high intensity" one ( Figure 2).
An extreme, but plausible, example of intensification is used for these comparisons. Specifically, an increase in colonies per apiary from 9 to 225 colonies, a change to a lattice configuration and a 10-fold increase in between-colony infection (0.015-0.15 per bee per day), demonstrated in Figure 2. The difference in the R 0 before and after intensification is how we estimate "additional R 0 ". This permits the interaction (nonadditive) effects of our three aspects of intensification. The "additional R 0 " can then be used in combination with the analytically derived relationship between R 0 and prevalence (see model (1) and Equations (2a-c)) to characterize how intensification affects disease prevalence. We focus on disease prevalence as both models show rapid pathogen spread across apiaries, such that infection prevalence at the endemic equilibrium was the major result differentiating modelling scenarios (S.I. Figures S4 and S5).

| RE SULTS
Our main results constitute three main characterizations of this system: the relationship between R 0 and pathogen prevalence; the effects of intensification on R 0 ; and by combination of these relationships, the effect of intensification on pathogen prevalence. The relationship between R 0 and pathogen prevalence is principally derived from the analytical model (presented first in these results) but is confirmed to broadly agree with the agent-based model (presented second). The relationship between intensification and R 0 is principally derived from the ABM, presented second, but is partly explored in the analytical model presented first. The critical overall result is the combination of these relationships, presented last and visualized in Figure 5, demonstrating how intensification impacts disease prevalence. Detailed derivation, exploration and testing of both models is detailed in the Supplementary Information.
Both model (1) and the ABM simulations show that, for a given number of colonies per apiary, R 0 is always greatest for the lattice arrangement-the most highly connected configuration. As the number of colonies per apiary increases (increasing n), the values of R 0 in both the array and lattice configurations increase (Figure 3a,b), while the R 0 for the circular configuration remains unchanged (see R 0 equations). The increase in R 0 from the addition of colonies asymptotes quickly due to convergence in the mean number of neighbours across the apiary; this is also why the R 0 for the circular apiary is Illustrative schematic of the "intensification" treatment as it is used in parts of this manuscript. We show the apiary used to estimate "base R 0 " (left) compared to the intensified apiary (right) reflecting an increase in number of colonies from 9 to 225, a change from an array to a lattice, and a 10-fold increase in movement of honeybees between colonies (illustrated using arrow weight) from a likelihood of 0.015 per bee per day to 0.15. Note that for the intensified apiary, not all 225 colonies are shown, with missing colonies denoted by ellipses (...) independent of number of colonies as the number of neighbours per colony remains two. This explains why R 0 for an array arrangement approaches the R 0 value for a circular arrangement as the number of colonies increases.
If R 0 > 1, the pathogen will rapidly invade (see S.I. Section 1, Figure S5) and each colony will reach a stable population size and infection prevalence, called the endemic equilibrium (See S.I. Section 1). Mathematically the disease prevalence at equilibrium for colony j is I j */(I j *+S j *), where S j * is the number of susceptible honeybees and I j * is the number of infectious honeybees in colony j at equilibrium.
The endemic equilibrium for the circular configuration model can be solved explicitly (see S.I. Section 1). Due to symmetry, all colonies within the circular apiary have disease prevalence at the endemic equilibrium of: We can approximate the endemic equilibrium for the lattice and array configured models using perturbation theory, assuming 0 < b « 1 (See S.I. Section 1). The approximate disease prevalence in colony j at equilibrium for a colony in the array or lattice configurations is: where l is the number of neighbours that colony j has. For any given set of parameters, we can, therefore, formulate both R 0 and prevalence, allowing us to characterize the relationship shown in Figure 3c.
We show analytically, and in the ABM (S.I. Section 3) that intensification in the form of an increase in colonies or an increase in movement between colonies increases R 0 (Figure 3a, b). Figure 4b shows the additional R 0 caused by our most extreme plausible changes in apiary management. The change in R 0 caused by increasing apiary size rapidly asymptotes (Figure 3a, b).
The effect of intensification is dependent on the base R 0 -for small base R 0 , intensification causes little additional R 0 , but at intermediate or high base R 0 , intensification leads to large additional R 0 ( Figure 4b). While the increase in R 0 is largest for an already large base R 0 , this relationship saturates and the relative increase in R 0 for a given base R 0 stays relatively constant for large base R 0 values. The relationship shows a strong nonlinearity when examining all three aspects of intensification in combination.
By understanding the effect of intensification on R 0 ( Figure 4b) and by characterizing the relationship between R 0 and Relationships between number of colonies, R 0 , and prevalence from model (1). Figure 3a and 3b demonstrate that the effect on R 0 for different degrees of intensification rapidly asymptotes, justifying our "single intensification" treatment ( Figure 2). Figure 3c defines the relationship between R 0 and prevalence, the shape of which critically determines our main result (see Figure 5). Technical description: (a) When R 0 = 30 for a single-colony apiary, the addition of colonies yields a maximum increase in R 0 of 12.7 for the lattice and 4.5 for the array. (b) When R 0 = 2 for a single colony, there is a maximum increase in R 0 of 0.85 for the lattice and 0.29 for the array, when colonies are added. Recall that the R 0 for the circle is independent of n (see (2b)), and hence absent from the figure. Parameter values are set to: v = 0.1, m = 0.0272, = 1,600 and in a) a + b = 4.32485 × 10 -6 and in b) a + b = 6.48725 × 10 -5 . The transmissibility is what affects base R 0 . Black dots are values where between-colony transmission is held at 10% of total transmission, with the bottom and top of the bars representing 1% and 20% of the total transmission being between hives, "b", respectively. (c) The relationship between R 0 and disease prevalence. The range of R 0 values is generated by varying the overall transmission rate (i.e. a + b) from 2.143 × 10 -6 to 1.178 × 10 -4 as reported by Roberts and Hughes (2015) for Nosema ceranae disease prevalence (Figures 3c and 3a), we can show how intensification impacts disease prevalences. We approximate the nonlinear relationship between "base R 0 " (pathogen phenotype) and the "additional R 0 " (effect of intensification) in Figure 4b. We use a bootstrapping approach to create 1,000 subsamples (subsample size = 10% of full sample with replacement) of our combined approach. Each subsample is used to generate a nonlinear model of the form y = ax/ (b + x c ), where y is "additional R 0 " and x is "base R 0 ", using a nonlinear least squares approach in R (v 3.3.1). The relationship generated using the full sample is plotted in Figure 4b.
We combine this relationship characterizing how base R 0 affects intensified additional R 0 (Figure 4b) with the derived relationship between R 0 and pathogen prevalence shown in Figure 3c, allowing us to predict how intensification impacts prevalences ( Figure 5). Figure 5a shows the proportion of bees infected by a given (base R 0 ) pathogen for the two apiaries in Figure 2. The difference in disease prevalence between these lines is the impact of intensification and is plotted in Figure 5b. Figure 5b shows a distinctly peaked relationship between base R 0 and the impact of intensification, with the impact of intensification peaking around base R 0 = 3.3, and then rapidly declining.
Even at its peak, the effect of intensification (which is as extreme as plausible), leads to an additional ~18% of bees infected at disease equilibrium. We present Figure 5 as the most important graphic for understanding the overall conclusions of this study, as the apparent "small" shift in R 0 required to double prevalence (Figures 3c and 4a) is actually very difficult to achieve for low R 0 pathogens (see Figures   3b and 4b), resulting in the "maximum plausible" change shown by the peak in Figure 5b (~18.5%).
We contextualize these results by calculating an estimate of the lower-bound of R 0 value for a honeybee pathogen (see highlighted regions in Figure 5). We identified this region based on empirical data for the microsporidian pathogen Nosema ceranae; this was the only pathogen for which experimentally derived transmission rates as well as robust information on mortality due to infection is avail- -Hernández et al., 2011;Paxton, Klee, Korpela, & Fries, 2007;Roberts & Hughes, 2015). To estimate the plausible R 0 boundary in our model for this pathogen, we parameterized our mathematical model using the lowest empirically supported transmission value with the highest supported additional mortality, and fixed movement of honeybees between colonies at its lowest supported natural rate (Currie & Jay, 1991). We then calculated the R 0 for a circular apiary due to its scale independence.

| D ISCUSS I ON
Our results present a counterintuitive picture of apicultural intensification and its consequences on disease prevalence within apiaries.
Even in their most plausibly extreme cases, changes in the number of colonies, their spatial arrangement and transmission rates between colonies (reflecting management intensification  had only a small effect on the severity of disease at the apiary level for pathogens of interest. Apicultural intensification leads to large gains in R 0 when R 0 is initially high and small gains in R 0 when R 0 is initially low (Figure 4b). However, increases in R 0 cause large increases in prevalence only when R 0 is initially low (Figures 3c and   4a). Pathogens with a base R 0 ≈ 3 benefit most from intensification in terms of increased prevalence ( Figure 5); As discussed below, we argue that there is likely to be a high base R 0 in important honeybee diseases and, therefore, our models suggest that there is likely to F I G U R E 4 Results from the ABM. Figure 4a demonstrates the agreement between the ABM and analytical model; Figure 4b presents the critical relationship estimated from the ABM relating base R 0 to the increase in R 0 following intensification (see Figure 2), the shape of which critically determines our main result (see Figure 5). Technical description: (a) shows agreement between the stochastic simulations (ABM) and analytical model ( Figure 3c); using the following equivalent model parameterization to that for Figure 3c: Circular configuration, n = 9, M = 58,200, Φ = 1,600, 5 × 10 -6 ≤ β ≤ 1 × 10 -4 , ν = 0.1, ρ = 0.1 (see S.I. Section 2). (b) Examines how an extreme example of intensification (see Figure 2) alters R 0 across a range of different "base R 0 " values determined by pathogen phenotype using the ABM. Grey points represent individual simulation comparisons, black points represent mean values. Base R 0 values are unevenly distributed across the range due to R 0 being an emergent property of the system in both plot panels. We derive a nonlinear relationship between "base R 0 " and "additional R 0 " for panel b, corresponding to Figure  However, if a pathogen emerges with a relatively low R 0 , our model does indicate that extreme intensification could lead to a significant increase in prevalence of approximately 18.5%. Therefore, if intensification increases the risk of novel pathogen emergence, then these newly emerged pathogens would benefit from intensification, as it would significantly increase their disease prevalence, relative to the pre-intensified apiary.
Our models most closely resemble the ecology of a directly transmitted microparasite able to infect individual honeybees at any life stage, conceptually similar to the microsporidian pathogens Nosema spp. (Fantham & Porter, 1912). Nosema is a major concern to beekeepers world-wide (Higes et al., 2008(Higes et al., , 2009Paxton, 2010), and has a minimum estimated base R 0 of 23 ( Figure 5) when modelled here. We found that apicultural intensification, in the context of a pathogen with an initial R 0 of 23, leads to a maximum 6.6% increase in disease prevalence. Our models predicted disease prevalences of up to 90% (Figure 3, Figure 5; S.I. Section 3), which while high, are empirically supported for the honeybee system (Higes et al., 2008;Kielmanowicz et al., 2015), and feature in other modelling studies that use similar transmission parameters to ours (Betti, Wahl, & Zamir, 2014). Nosema was the only pathogen for which there are direct empirical studies characterizing its transmissibility, however, other honeybee pathogens such as deformed wing virus are also well studied. While estimating an R 0 for DWV is difficult due to active management by beekeepers, maximum reported prevalences that may be indicative of its true "unmanaged" R 0 are high, for example 73% in Natsopoulou et al. (2017), 80% in Budge et al. (2015) and 100% in Stamets et al. (2018). These high prevalences are consistent with high R 0 values (Figures 3c and 4a, and S.I. (Section 3)).
We additionally explored the behaviour of a more specific model, using an age-structured approach to infection dynamics, where only larvae are vulnerable to infection and develop into infectious adults with a high pathogen-associated mortality (as might be appropriate for pathogens such as the acute paralysis virus complex (Martin, 2001)), presented in the S.I. (Section 3). Convergence to equilibrium happens more slowly than the main model presented here, but still occurs quickly (within a single beekeeping season; see S.I. 3 Figure   S7). However adult-bee infection prevalence is far lower than seen in our SI model (S.I. Figure S7)-this is in agreement with observations of lower prevalence of paralysis viruses (Budge et al., 2015). Notably, the endemic equilibrium prevalence increases only by small magnitudes as movement between colonies or apiary sizes are drastically increased (S.I. Figure S7), in agreement with our main general result.
This equivalence in behaviour between different models reflecting large disparities in infection mechanics and different endemic prevalences demonstrates that these results are likely generalizable to many honeybee pathogens.
We find rapid spread of a given pathogen across an apiary, which quickly reaches endemic equilibrium (S.I. Figures S4-S6). While pathogens with a higher R 0 reach this equilibrium more quickly, there is universally rapid spread. Given this result, we mainly focussed on the disease prevalence experienced at endemic equilibrium. Despite assuming transmission only to nearest neighbours, pathogen spread occurs rapidly, and the nearest neighbour assumption alters this very little when removed or relaxed (see S.I. Figure S6). The rate at which epidemics are established in our model is also in agreement with other honeybee pathogen models. For example, Jatulan, Rabajante, Banaay, Fajardo, and Jose (2015) show that a single infectious adult causes an American Foulbrood (Paenibacillus larvae) epidemic that peaks within 50 days. Whilst they do not explicitly find an R 0 for P.
larvae, the short timescales characterizing their epidemics are in line with ours (S.I. Section 3), suggesting high R 0 values and that their model would behave similarly to ours at an apiary scale.
F I G U R E 5 Depictions of our critical finding characterizing the maximum (peak), and likely (shaded region), increases in prevalence of a pathogen following local intensification of apiculture. High prevalence even in "low intensity" (see Figure 2) systems yields little opportunity for large increases in prevalence. Panel (a) shows the proportion of bees infected (prevalence) in non-intensified apiaries (lower red line) compared to intensified apiaries (upper blue line), take from the mean values derived in Figure 4b and the relationship shown in Figure 3c. The shaded grey area between these curves is the additional prevalence caused by intensificationthe "impact of intensification". This is plotted in panel (b) where the black line represents the mean relationship, and the grey lines represent 1,000 bootstrapped samples. The vertical dashed line and yellow-shaded region of the graphs to the right of the dashed line show a lowest estimated value of R 0 for Nosema ceranae. Figures start at R 0 = 1.0008 Our intercolony transmission can be understood to capture multiple processes arriving from beekeeper management such as brood transplantation or reduced distance between colonies  as well as recognized transmission routes such as honeybee drift (Jay, 1965). Our approach was informed by studies which have focussed on how changes in the number of colonies and apiary configurations (Jay, 1966(Jay, , 1968 alter drift (Dynes et al., 2017). Links between drift-mediated pathogen transmission and colony numbers have been documented for a variety of pathogens (Seeley & Smith, 2015)-including brood specialized and non-specialized, micro-and macro-parasites (Belloy et al., 2007;Budge et al., 2010;Dynes et al., 2017;Nolan & Delaplane, 2017). Larger numbers of colonies per apiary are a driver of higher drift (Currie & Jay, 1991), as are changes in apiary arrangement (Dynes et al., 2019;Jay, 1966). While beekeepers typically maintain equal distances between their colonies regardless of how many colonies are in the apiary (such that larger apiaries have a bigger area footprint), our approach of increasing between-colony transmission in larger apiaries would also capture any additional transmission from spatial crowding.
Two clear candidates for future development of this model include seasonality and demography, which are closely linked.
Honeybee demography within a colony influences epidemiology (Betti, Wahl, & Zamir, 2016) due in part to the temporal polyethism of task allocation influencing exposure and immunity (Calderone & Page, 1996), as well as the flexible ability of honeybees to regain immune function when they revert roles (Amdam et al., 2005;Robinson, Page, Strambi, & Strambi, 1992 (Fürst et al., 2014;Manley et al., 2015;McMahon et al., 2015McMahon et al., , 2018. Honeybee colony density across a landscape, therefore, has implications for wild pollinator health (Cohen et al., 2017;Graystock et al., 2016), however, our results suggest that increased stocking of honeybees may have smaller impacts on local pollinator infectious disease dynamics than may have been previously thought.
Other industrialized agricultural livestock systems reflect extreme host densities similar to those in this study. However, the R 0 for honeybee diseases may exceed that of other livestock diseases.
We compare our lower threshold estimate for the R 0 of N. ceranae to all available R 0 values for livestock diseases that we could readily find in the literature ( Figure S9, see S.I. Section 4). Notably, all other livestock diseases for which R 0 estimates exist show minimum R 0 values far below our honeybee estimate, however, examples of agricultural R 0 values as high or higher than those we present for honeybees do also exist. There is, therefore, a clear need to develop explicit models of agricultural intensification scenarios for important agricultural disease.
Overall, our findings represent the first stage in developing robust epidemiological models for studying honeybee pathogens at an apiary scale. In the face of increasing challenges to global apiculture, our models predict that the size of apiaries per se is not causing no- Finally, this study demonstrates that conventional thought on how agricultural intensification influences disease may not be robust in the face of system-specific ecological nuance.

DATA AVA I L A B I L I T Y S TAT E M E N T
The agent-based model is made available in association with this manuscript via Dryad Digital Repository https ://doi.org/10.5061/ dryad.rn2j5p0 .