The "Brownian bees" model describes a system of N-independent branching Brownian particles. At each branching event the particle farthest from the origin is removed so that the number of particles remains constant at all times. Berestycki et al. [arXiv:2006.06486] proved that at N→∞ the coarse-grained spatial density of this particle system lives in a spherically symmetric domain and is described by the solution of a free boundary problem for a deterministic reaction-diffusion equation. Furthermore, they showed [arXiv:2005.09384] that, at long times, this solution approaches a unique spherically symmetric steady state with compact support: a sphere whose radius ℓ_{0} depends on the spatial dimension d. Here we study fluctuations in this system in the limit of large N due to the stochastic character of the branching Brownian motion, and we focus on persistent fluctuations of the swarm size. We evaluate the probability density P(ℓ,N,T) that the maximum distance of a particle from the origin remains smaller than a specified value ℓ<ℓ_{0} or larger than a specified value ℓ>ℓ_{0} on a time interval 0<t<T, where T is very large. We argue that P(ℓ,N,T) exhibits the large-deviation form -lnP≃NTR_{d}(ℓ). For all d's we obtain asymptotics of the rate function R_{d}(ℓ) in the regimes ℓ≪ℓ_{0},ℓ≫ℓ_{0}, and |ℓ-ℓ_{0}|≪ℓ_{0}. For d=1 the whole rate function can be calculated analytically. We obtain these results by determining the optimal (most probable) density profile of the swarm, conditioned on the specified ℓ and by arguing that this density profile is spherically symmetric with its center at the origin.