Radio-Frequency Electromagnetic Field Exposure of Western Honey Bees

Radio-frequency electromagnetic fields (RF-EMFs) can be absorbed in all living organisms, including Western Honey Bees (Apis Mellifera). This is an ecologically and economically important global insect species that is continuously exposed to environmental RF-EMFs. This exposure is studied numerically and experimentally in this manuscript. To this aim, numerical simulations using honey bee models, obtained using micro-CT scanning, were implemented to determine RF absorbed power as a function of frequency in the 0.6 to 120 GHz range. Five different models of honey bees were obtained and simulated: two workers, a drone, a larva, and a queen. The simulations were combined with in-situ measurements of environmental RF-EMF exposure near beehives in Belgium in order to estimate realistic exposure and absorbed power values for honey bees. Our analysis shows that a relatively small shift of 10% of environmental incident power density from frequencies below 3 GHz to higher frequencies will lead to a relative increase in absorbed power of a factor higher than 3.

information on the evolution of such absorption as a honey bee goes through different developmental stages. Nor is it clear whether this RF absorption is realistic for other castes, such as drones or queens, in a bee colony.
Therefore, the aims of this study were to numerically evaluate RF-EMF absorption in western honey bees and validate the frequency dependency of this absorption during various developmental stages and experimentally quantify real-life exposure of bees. To this aim, numerical simulations were executed to determine the absorption of RF-EMFs in five different honey bee models: a larva, a queen, two workers, and one drone, obtained using micro-CT imaging. These simulations were implemented as a function of frequency in a broad band, 0.6 GHz up to 120 GHz, that can be used to model both current and future telecommunication frequencies. In parallel, RF-EMF exposure measurements were executed near five bee hives in Belgium, in order to quantify the real exposure of such honey bees. Finally, these measured values were used to rescale the numerical simulations in order to quantify real honey bee absorption and assess a potential change in absorption in case a shift in operation frequencies in future telecommunication networks would occur.

Methods
Studied honey bees, imaging technique, and model development. Images of the studied insects are shown in Fig. 1. All studied insects are western honey bees (Apis mellifera), which is the most commonly used honey bee worldwide. Honey bees within a colony are subdivided into different castes. An active viable honeybee colony contains only one queen bee who spends most of her time laying 2,000 to 3,000 eggs per day. The queen is the only reproductive female within the colony and her health is vitally important to the survival of her colony. Damage to her ovaries has the potential to effect the function and survival of her progeny. A queen typically lives between approximately three and five years. From early spring time to mid-summer the queen lays unfertilized "haploid" eggs which develop into drone bees. All drones are males. Their specific role is to mate with a virgin queen so that she can initiate the propagation of a new colony. During this mating season, there are approximately 3,000 to 5,000 drones within any given colony. Drones typically live between one to two months.
A healthy honey bee colony can contain approximately 50,000 individuals. Most of these are sterile, female, worker bees. Worker bees perform all the tasks within a colony to keep it full of provisions and free from disease. This involves feeding and nursing larvae, foraging for nectar and pollen, storing nectar and pollen, guarding the entrance, tending to the hygiene of the queen-workers-drones and maintaining a clean hive environment. Workers live for three to four weeks during the active seasons (spring-summer-autumn) and approximately three months during the colder inactive season (winter). There are approximately 3,000 (winter) to 10,000 (summer) larvae present at any given time.
We chose representatives from all three castes within a honeybee colony, one queen bee, two worker bees, one drone bee and one worker larva. All honey bees were scanned at the Western Sydney University National Imaging Facility (Sydney, Australia) using a bench-top MicroCT scanner (Quantum GX MicroCT Imaging System, PerkinElmer, Hopkinton, MA, USA). The parameters used during this scanning depended on the scanned bee. Such scans are made using different projections, at different time intervals on the scanners settings. The rotation between projections also depends on the scanner's settings and the studied honey bee (see below for full description).
Worker 1. The insect named 'Worker 1' is the same bee studied in 11 , which had a full body length of approximately 11.0 mm long, is 5.0 mm wide, and had a mass of approximately 900 mg. During the scanning of Worker 1, the Micro-CT scanner was operated using the following parameters: 50 kVp, 80 mA, and a 2048 × 2048 pixels image matrix. This resulted in scans with a 20 μm isotropic voxel size. Each projection had a scanning time of 3.0 s, with 3.0 s rotation time in between projections. The total scan time for Worker 1 was approximately 18 min.
Worker 2. The second honey bee worker (Worker 2) has a full body length of 13 mm with cross sectional dimensions of 6.8 mm and 5.4 mm and a mass of approximately 900 mg. For Worker 2, the scanner was operated using the following parameters: 40 kVp, 70 mA, and a 2048 × 2048 pixels image matrix. The isotropic voxel size was 100 μm. Each projection had a scanning time of 1.5 s. There was a 3.0 s rotation time in between each projection. The total scan time for the whole bee was approximately 10 min.
Larva. Larvae of this age (three weeks) are typically approximately 16 mm long with an approximate mass of 900 mg. The scanned larva was curled up, which made estimating its full body dimensions difficult, but the sample fitted within a 14 × 7 × 15 mm 3 box. This scanning of the larva was done using the following parameters: 50 kVp, 80 mA, and a 2048 × 2048 pixels image matrix. This resulted in scans with a 20 μm isotropic voxel size. Each projection had a scanning time of 3.0 s. and with a 3.0 s rotation time this resulted in a total scan time for the larva of 18 min.
Male drone. The drone has a full body length of 18 mm with cross sectional dimensions of 7.2 mm and 9.4 mm and an approximate mass of 1 g. During the scanning of the drone, the Micro-CT scanner was operated using the following parameters: 40 kVp, 70 mA, and a 2048 × 2048 pixels image matrix. The isotropic voxel size was 100 μm. Each projection had a scanning time of 1.5 s. The full scan took 180 projections and there was a 3.0 s rotation time in between each projection. The total scan time for the whole bee was approximately 10 min.
Queen bee. The QB has a full body length of 19 mm and cross sectional dimensions of 7.5 times 7.1 mm 2 and an approximate mass of 1100 mg. The queen was scanned was using the following parameters: 40 kVp, 70 mA, and a 2048 × 2048 pixels image matrix. The isotropic voxel size was 250 μm. Each projection had a scanning time of 1.5 s. There was a 1.5 s rotation time in between each projection. The total scan time for the queen bee was approximately 10 min.
Development of 3D models. The software running on the Quantum GX, bench-top MicroCT scanner was used for all honey bees to reconstruct the 180 projection images. Those were then converted into a 2D rendered image stack of 512, 16 bit bitmap images. Finally, the BeeView volume rendering software (DISECT Systems Ltd, Suffolk, UK) was used to acquire Bee volume data from the image stack. All 3D models of the insects were created using the software TomoMask (www.tomomask.com). We used the same approach as in 11 . The image stack for each honey bee was imported into TomoMask, which also required the pixel and slice spacing. The software generated a 3D model using a marching cubes algorithm 21 . This model was then exported as an STL (STereo Lithography) 22 file. This is a commonly used format to describe surface geometry. The models were also smoothed using the Taubin λ/μ smoothing scheme 23 implemented in MeshLab 24 . The dimensions of the models and mesh integrity were checked (and corrected if necessary) before simulations using Netfabb (Autodesk, San Rafael, CA, USA).
Numerical simulations and RF EMF exposure conditions. Electromagnetic, numerical simulations were executed to estimate electromagnetic fields in and around the honey bees under far-field exposure. Far-field exposure is in this manuscript defined as RF-EMF sources being more than 2D 2 /λ away from the insects, with D the largest dimension of the RF source and λ the wavelength of the RF-EMFs. This is often referred to as the Fraunhofer far-field limit 25 . In general, far-field RF-EMF sources can be located in any direction from the honey bees. Therefore, different approaches exist to model such far-field exposure to RF-EMFs: a stochastic method where far-field exposure is decomposed in sets of plane waves according to certain statistics is used in 26,27 , while a more limited set of plane-wave exposures coming from six predefined directions along the main axis of the exposed subject or animal are considered in 11,28 . In this study, we have chosen to work with the latter method. We have modeled exposure of the studied honey bees by a set of 12 incident plane waves traveling along six directions defined by a Cartesian coordinate system, see Fig. 2. For each direction, two orthogonal incident electric field polarizations were chosen, since any other free-space E-field polarization can be obtained using a linear combination of both. All incident plane waves have a root-mean squared electric field strength of 1 V/m. This value is chosen to facilitate renormalization to any potential value of incident field strenght. Numerical simulations were executed using the Finite-Difference Time-Domain (FDTD) method implemented in Sim4life (ZMT, Zurich, Switzerland). This is a common technique used to determine RF-EMF in and near homogeneous and heterogeneous dielectric objects 11,26,28 , such as the honey bees studied in this paper. In this method, the simulation domain is divided in cubes using a three-dimensional rectilinear grid. Depending on the wavelength, feature sizes of the objects in the simulations, and the desired spatial accuracy, a different spatial step is used to discretize the simulation. The FDTD algorithm requires a grid step smaller than one tenth of the smallest wavelength in the simulation domain in order to return stable solutions 29 . Since this is a time-domain technique, it requires a predefined simulation time in order to reach a steady-state solution, which will again depend on the chosen spatial resolution, the wavelength, and the size of the simulation domain.
www.nature.com/scientificreports www.nature.com/scientificreports/ We executed numerical simulations at nine harmonic frequencies from 0.6-120 GHz (sinusoidal waves at a single frequency). The lower and upper frequency limits were chosen because they correspond to the current limits in terms of simulation size and length that can realistically be supported by our simulation hardware. The simulated frequencies are listed in Table 1 alongside the chosen grid steps in the simulation domain and the number of periods used for every simulation. These settings were the same for each of the five studied honey bee models. The studied insects have certain dielectric properties, quantified using the relative permittivity (ε r ) and conductivity (σ). We did not measure the dielectric properties of the studied insects. Instead, we assigned dielectric parameters obtained from 11 . The value at 1 GHz is obtained using the same literature database and interpolation presented in 11 . Table 1 lists these properties. All insects were modeled as homogeneous objects. These configurations resulted in 12 (plane waves) × 9 (frequencies) × 5 (honey bees) = 540 simulation results.
After each simulation, the internal electric field in the insect model was extracted and used to calculate the total absorbed RF-EMF power (P abs ) in the honey bee. P abs is calculated as the integrated product of the conductivity and the squared internal electric field strenght (E int ) over the total volume (V) of the insect: We report P abs rather than specific absorption rate (SAR) values since we did not measure the mass and density of all the simulated honey bees. P abs is an important quantity since dielectric heating of an insect is proportional to absorbed RF-EMF power 2 .
In order to validate our simulations we tested the influence of four simulation settings on the RF-EMF P abs : grid step size, dielectric parameters, angle of incidence, and number of simulated periods. The influence of the grid step is expected to be the most significant at the highest simulated frequency (120 GHz), since the chosen

GHz 1.2 GHz 2 GHz 3 GHz 6 GHz 12 GHz 24 GHz 60 GHz 120 GHz
Maximal grid step (mm) www.nature.com/scientificreports www.nature.com/scientificreports/ maximal grid step of 0.05 mm is closest to the smallest wavelength in the simulation domain at that frequency in the tissue (0.05 mm = 0.045 λ). Therefore the maximal grid step was set to 25 μm for exposure configuration number 2 in Fig. 2 for both the Larva and Worker 2 phantoms. In 11 , it was demonstrated that the maximal uncertainty on the dielectric parameters occurs between 2 and 3 GHz, with maximal relative deviations of 40%. In order to test the dependency of our simulation results on the chosen dielectric parameters, we executed four additional FDTD simulations in exposure configuration number 2 shown in Fig. 2 using the Worker 2 phantom. In these simulations the dielectric parameters (ε,σ) were changed to: (1.5.ε, 1.5.σ), (0.5.ε, 1.5.σ), (1.5.ε, 0.5.σ), and (0.5.ε, 0.5.σ), respectively, allowing for a potential 50% deviation on the dielectric parameters, which should be larger than the uncertainty on the chosen dielectric parameters. We chose to model RF-EMF exposure of the studied honey bees using plane waves incident from 6 directions. However, it is uncertain whether this set of plane waves provides a complete overview of the full range in P abs as function of the angle of incidence. In order to validate our exposure set up, we have executed 20 additional FDTD simulations at 6 GHz using the Worker 2 phantom, where the elevation, azimuth, and polarization angles were generated according to uniform distributions between [0, π], [0, 2π], and [0, 2π], respectively. The settings of these FDTD simulations were the same as those shown in Table 1. Finally, the number of simulated periods was tested at 120 GHz for the Worker 2 phantom in exposure configuration number 2 shown in Fig. 2 by increasing the number of simulated periods to 120 instead of 30, see Table 1. After each of these validation simulations, the P abs was extracted and compared to the one obtained in the original simulation set.

RF-EMF field measurements.
In order to quantify current RF-EMF exposure of honey bees in real exposure scenarios, we executed RF-EMF exposure measurements at five sets of bee hives in Belgium at: Aalter, Merelbeke, Eeklo, Zomergem, and Drongen, see Fig. 3(a). At each measurement site, three different measurements were executed in order to quantify RF-EMF exposure.
First, a spectrum analyzer of the type FSL6 (R&S Belgium, Excelsiorlaan 31 1930 Zaventem Belgium) connected to a triaxial isotropic antenna was used to perform a broad-band RF overview measurement from 80 MHz to 6 GHz. These measurements were executed in two steps: first spectral overview measurements were executed from 0.08-3 GHz using a tri-axial antenna TS-EMF (Rhode and Schwartz, dynamic range of 1 mV/m-100 V/m for the frequency range of 80 MHz-3 GHz), followed by measurements from 3-6 GHz using a Clampco AT6000 antenna. At one out of five measurement sites, Drongen, a conical dipole antenna PCD 8250 (Seibersdorf Laboratories, Seibersdorf, Austria) was used for the 80 MHz -3 GHz measurements. This antenna was rotated to obtain three orthogonal polarizations of the electric field. During these overview measurements, the spectrum analyzer measured in maximum-hold modus during 17 and 9 minutes in the lower and higher frequency bands, respectively. The antennas were supported by a plastic tripod and were placed at 1 m in front of the bee hive at a height of 1.5 m from the ground level. Figure 3 shows the studied bee hives and the measurement set up in the field. The 1.5 m height is a typical height at which such EM field measurements 30 . Additionally, this height is mentioned in the ECC(02)04 standard 31 . The purpose of these measurements was to get an overview of which frequency bands were in use at the respective sites. These frequency bands were then investigated further in the second measurements.
Second, the same spectrum analyzer was connected to the tri-axial antenna TS-EMF which was again supported by the same tripod at a height of 1.5 m. The tripod was placed at two distances of 1 and 2 m from the central bee hive. The spectrum analyzer performed root-mean square electric field strength (E RMS ) measurements over a measurement period of 6 minutes 2 in each of the telecommunication frequency bands identified using the first measurement. Each of the three electric field components (E x , E y , E z ) were measured individually. E RMS was then obtained as the square root of the sum of squares of the individual components.
The spectrum analyzer measurements in terms of received power on the antenna were then recalculated using the known antenna factor of the tri-axial antenna to incident root-mean-squared electric field strength. The E RMS,i values in each frequency band (i) were then summed quadratically and the square root of that sum is listed as the total instantaneous electric field strength (E RMS,tot ).
The measurement procedure and measurement settings for these RF-EMF exposure measurements are presented in 32 . The expanded measurement uncertainty (95% confidence interval) for electric field strength measurements using this set up is ±3 dB 30 .This measurement setup enables the most accurate assessment of in situ exposure from various RF-EMF sources 30 .
Third, a broadband exposure measurement was executed using a Narda NBM-550 probe (Narda, Hauppauge,NY, USA) connected to an EF 0691 broad-band probe (Narda, Hauppauge, NY, USA) which has a frequency span from 100 kHz to 6 GHz, thus including so-called intermediate frequencies (IF). These IF fields are not considered in our numerical simulations. However, we measured those to provide a complete overview of the exposure to electromagnetic field below 6 GHz. The NMB probe was placed on top of the central bee hive and was left there during both RF measurements. The device measured and registered root-mean-squared electric field strengths with a period of 1 s. From those time series of measurements, we obtained the time average and the maximal value. (2020) 10:461 | https://doi.org/10.1038/s41598-019-56948-0 www.nature.com/scientificreports www.nature.com/scientificreports/ The researchers that executed the RF-EMF field measurements did not use personal devices during the measurements. All wireless devices brought to the measurement site by the researchers were operated in flight mode, i.e. any wireless transmissions by those devices were not allowed.

Estimation of realistic RF-EMF absorbed power in honey bees.
Realistic P abs absorbed in honey bees can be obtained by rescaling the simulated P abs values using the measured incident field strengths. Therefore, we linearly averaged the total E RMS values measured near the five bee hives at two different positions to obtain an average E RMS,avg value. In order to estimate exposure of honey bees in current wireless networks, we averaged the P abs values using: with f i = 0.6, 1.2, 2, 3 GHz. We only considered P abs values < 3 GHz, since our measurements will show that there are only incident RF-EMFs below 3 GHz in the current environment of honey bees in Belgium. This value is then rescaled using:  www.nature.com/scientificreports www.nature.com/scientificreports/ In order to estimate the effect of a fraction (p ∈ [0, 1]) of the RF-EMF incident fields shifting to frequencies higher than 3 GHz we also determine the average P abs for frequencies higher than 3 GHz, using: abs av j abs j , 1 with f j = 6, 12, 24, 60,120 GHz. The realistic P abs,real (p) for a fraction p of the power shifted to frequencies higher than 3 GHz is then calculated as:

Results
Numerical simulations. Figure   www.nature.com/scientificreports www.nature.com/scientificreports/ field strengths decreases very rapidly within the phantom and electric fields are basically only present in the outer layers of the insect. This is caused by a decrease in skin depth that is driven by the increase in conductivity at higher frequencies, see Table 1. Note that the total RF-EMF absorbed power in the insect scales both with the internal electric field strength and the conductivity. Figure 5 shows the normalized RF-EMF P abs as a function of frequency for the five studied insects from 0.6 GHz up to 120 GHz. The curves connect the linear averages of the 12 P abs values obtained for each honey bee at each simulated frequency, while the whiskers indicate the minimum and maximum P abs values found at those frequencies. All P abs values are normalized to an incident field strength of 1 V/m. Figure 5 shows an increase of P abs over frequency for all studied phantoms up to 6 GHz. When comparing the average P abs at 0.6 GHz and 6 GHz, we found relative increases of factors of 16, 35, 72, 121, and 54 for the Worker Bee 1, Worker Bee 2, Drone, Larva, and queen Bee, respectively. The P abs slightly decreases over frequency beyond 12 GHz for all the studied honey bees. When comparing P abs at 12 GHz and 120 GHz, we found relative decreases of 26%, 34%, 33%, 32%, and 34% for the Worker Bee 1, Worker Bee 2, Drone, Larva, and Queen Bee, respectively. The spread on the P abs values obtained at each individual frequency reduces from up to a factor of 13 below 12 GHz to smaller than a factor 2.5 beyond 12 GHz. Figure 5 shows a general increase of P abs with increasing volume and surface area of the studied insects. Previous studies on whole-body averaged absorbed RF power and specific absorption rate of humans have shown a dependency of these quantities on the absorption cross section, a quantity that scales with volume and/or surface area of an exposed subject. When the diagonals of the smallest rectangular brick that contain the insect phantoms are considered, the honey bee with the smallest diagonal, Worker Bee 1 with a diagonal of 13 mm has the overall lowest average P abs . The Larva, Queen Bee, and Drone all have associated diagonals of 22 mm and have similar average P abs values as function of frequency. The Worker Bee 2 has a diagonal that falls in between Worker 1 and the other insects of 16 mm and also has an average P abs that falls in between the curve for the smaller worker and the other honey bee models, see Fig. 5. We attribute he differences between the two Worker Bee phantoms mainly to the difference in size of both phantoms. The larger Worker Bee 2 phantom has a larger diagonal, surface area, and volume. This leads to a higher absorption cross section 33 and higher P abs .
The maximal P abs for the five studied insect models occurs at those wavelengths that are close to the double of this diagonal, which suggests an absorption peak around half a wavelength. The maximum P abs for the Larva model lies in between 3 and 12 GHz, i.e. in between 25 and 100 mm in terms of λ, while the diagonal of said bounding box is 22 mm for the phantom. For the other studied insect models the maximum P abs lies in between 6 and 24 GHz, i.e. in between 23 and 50 mm in terms of λ, with associated phantom diagonals ranging from 16 mm to 22 mm.
As mentioned in the Methods section, the influence of dielectric parameters was studied with simulations using Worker 2 at 2 GHz with altered dielectric parameters. These resulted in P abs values of 6.3 × 10 −10 W, 6.3 × 10 −9 W, 3.1 × 10 −9 W, and 1.8 × 10 −9 W, in comparison to 2.0 × 10 −9 W for an incident field strength of 1 V/m. This corresponds to relative deviations of −69%, +210%, +50%, and −10%. These deviations are significant but smaller than the full range of a factor of 5 we observed for the larva at 2 GHz as a function of changing incident angle and polarization. These relative differences are small in comparison to the differences we observe over frequency for the same phantom: a factor of 121 over frequency from 0.6 to 6 GHz.
At 120 GHz we find a deviation on P abs smaller than 0.1% when 120 simulation periods are executed in comparison to 30 simulation periods in configuration number 2 shown in Fig. 2 for the Worker 2 phantom. Indicating that the number of simulated periods is sufficient for these simulations. At the same frequency and in the same simulation configuration, a reduction of the grid step with a factor of 2 resulted in a P abs of 8.6 × 10 −8 W and 3.1 × 10 −7 W for the Worker 2 and Larva phantoms, respectively, while the regular simulations with 0.1 mm www.nature.com/scientificreports www.nature.com/scientificreports/ and 0.05 mm grid steps, respectively, resulted in P abs values of 8.4 × 10 −8 W and 3.1 × 10 −7 W for an incident field strength of 1 V/m. This corresponds to relative deviations of 0.3% and 0.5% for the Worker 2 and the Larva phantoms, respectively, indicating that the chosen grid step was small enough to result in stable numerical results.
The set of 20 incident plane waves with randomized angles of incidence and polarization at 6 GHz using the Worker 2 phantom resulted in an average P abs of 4.5 × 10 −8 ± 1.6 × 10 −8 W for an incident field strength of 1 V/m, while the set of 12 incident plane waves used to model far-field exposure results in an average P abs of 6.5 × 10 −8 ± 5.3 × 10 −8 W at the same frequency. The value are fairly close, which indicates that the set of 12 incident plane waves along the main axes is a good proxy for average exposure under a randomized angle of incidence and polarization. The set of twelve plane waves does seem to overestimate exposure at the higher percentiles, since they are significantly higher than those obtained using the random set of plane waves. Figure 6 shows an example of an RF-EMF overview measurement at one of the five studied bee hives (Aalter). Figure 6 shows the relative electric field strength, normalized to the maximally measured electric field strength. The different peaks correspond to several individual frequency bands that are used for telecommunication and broadcasting signals. These frequency bands were then measured individually using the same set-up with triaxial antenna and spectrum analyzer at two positions relative to the bee hive on each measurement site using the measurement procedure described in 32 . Table 2 lists the measured E RMS values at the five studied bee hives shown in Fig. 3. As all these measurement sites were rural, private areas, there were no uplink (emissions from a user device to the network) transmissions found. Downlink (DL, this is network to user communication) signals were found at all measurement sites. These signals were generated by three different mobile telecommunications providers in fourteen different frequency bands.  www.nature.com/scientificreports www.nature.com/scientificreports/ area. WiFi was not measured at three out of five measurement sites. Additionally, at all measurement sites, RF EMFs emitted by a pulsed radar or other wireless technologies used in aeronautical surveillance were observed. The E RMS value of RF EMFs emitted by a radar cannot be accurately measured without having the specifications of the radar. Therefore, we can only measure the peak value over the 6 min measurement interval. These fields were the highest in Merelbeke, where at position 1 peak E-field values of 0.017 V/m and 2.2 V/m were measured at 1.09 GHz and 1.3 GHz, respectively, while at position 2 peak E-field values of 0.02 V/m and 2.9 V/m were measured at at 1.09 GHz and 1.3 GHz, respectively.

RF-EMF field measurements.
In order to provide the readers with a complete overview of the exposure to EMF fields below 6 GHz at the chosen measurement sites, Table 3 lists measured values in the 100 kHz to 6 GHz range using a broadband field   www.nature.com/scientificreports www.nature.com/scientificreports/ probe. All the average values are higher than what is obtained from the frequency-selective measurements presented in Table 2, as should be the case since a broader band is considered. Table 2, one can rescale the P abs values shown in Fig. 5 in order to obtain a realistic estimate of the absorbed RF-EMF power in honey bees P abs,real . The third to eight columns of the top row of Table 4 list P abs,real assuming that all incident E rms = 0.06 V/m is uniformly distributed over the simulated P abs values lower than 3 GHz. These values range from 0.1 nW for Worker 1 until 0.7 nW for the Larva and Queen Bee. In each subsequent row, 10% of the incident power density is transferred to frequencies higher than 3 GHz. This causes an increase in the estimated P abs,real (p). In order to quantify this increase, the five columns to the right show the relative increase in P abs,real (p) as p increases from 0 to 1. A full shift of all RF-EMF power to frequencies higher than 3 GHz -without changing the incident field strength -would result in relative increases in absorbed power between a factors 24-48 for the studied honey bee models. Even a relatively small shift of 10% of the incident power density to higher frequencies will lead to a relative increase in P abs of a factor higher than 3, see Table 4.

Discussion
This study investigates RF-EMF absorption in Western Honey Bees as a function of frequency in the 0.6 to 120 GHz range. To this aim, we used five different models of different honey bees: two workers, a drone, a larva, and a queen. These models were obtained using micro-CT imaging and used for FDTD simulations. These were used to evaluate far-field exposure of honey bees. This far-field exposure is modeled as a set of plane waves at harmonic frequencies between 0.6 and 120 GHz. The numerical simulations resulted in P abs as a function of frequency for the different studied honey bees. These simulations were combined with real RF-EMF exposure measurements near bee hives in Belgium in order to estimate realistic exposure values for honey bees.
Micro-CT imaging is a technique that has previously been shown to accurately scan insects 35,36 . The models used in this study have resolutions between 0.02 mm and 0.25 mm, which is larger than the resolution of the micro-CT models using in 11 . Since the smallest grid step used in our simulations is 0.05 mm, the ideal resolution of the insect models would be smaller than that. The larger resolution of the scanning is not a problem for the stability of the FDTD algorithm, but more spatial resolution could be obtained with the same simulation settings. It is expected that the micro-CT models used in this study lead to a better estimation of P abs and the spatial distribution of the electric fields than approximate models such as ellipsoids or cylinders 37 .
The results of our numerical simulations, see Fig. 5, show an increase of P abs with frequency up to 6-12 GHz. Figure 4 illustrates the mechanism behind this increase: as the frequency increases the EMFs are less likely to diffract around the honey bees, that are relatively small in comparison to the wavelengths <6 GHz, and can penetrate further in the models, generating higher internal electric fields and consequently higher P abs values. Figure 4 also shows why the whole-body averaged P abs does not increase beyond 12 GHz. As the conductivity increases, see Table 1, the electric fields will decay faster within the honey-bee phantoms, which leads to larger relative volumes within the insect with lower fields, see Fig. 4, which will also contribute to the whole-body averaged P abs . This effect also causes the P abs to have a smaller dependency (variation) on incident angle and polarization, see Fig. 5. We also observe that both the frequency-dependency of the P abs , i.e. the transition point between sharp increase in P abs over frequency and slight decrease over frequency, and the magnitude of the P abs , i.e. the offset of the P abs curve, depend on the honey bee's size. This effect was previously observed in 11 . In general, the results presented in this manuscript are in excellent agreement with those presented in 11 . The results in terms of P abs obtained for the honey bees in this study fall right in between those obtained in 11 for the smaller Australian Stingless Bee and the larger Desert Locust, which confirms again the dependency of P abs on phantom size. The same size-related effect was described for humans in 28,33,38 and comparable frequency trends were observed in humans that have larger full-body sizes at MHz frequencies 28,38 . It should be noted that this manuscript focused on exposure of individual insects in free space. In reality, honey bees might cluster, creating a larger absorption cross section and potentially higher absorption at lower frequencies.   www.nature.com/scientificreports www.nature.com/scientificreports/ The FDTD simulations presented in this manuscript use dielectric properties that were obtained from the literature survey executed in 11 . Ideally, these dielectric parameters would be obtained for the honey bees studied in this manuscript. However, as shown in 11 , most studies on dielectric properties of insects in literature 3,[39][40][41] show similar frequency dependencies of those dielectric parameters. We have executed additional numerical simulations to test for the uncertainty on the dielectric parameters and found deviations up to 210% on P abs , which is significant but still smaller than the variations that exist due to changing angle of incidence and polarization at a fixed frequency, or changes in frequency. We modeled the insects as homogeneous dielectric objects, while in reality they have heterogeneous dielectric parameters. Even though the FDTD algorithm will always require an averaging of dielectric parameters over the cube size, further developments in honey bee and insect phantoms should be focused on the inclusion of multiple tissues in order to refine these models.
In-situ RF-EMF measurements were executed using a measurement set up consisting out of a spectrum analyzer connected to an isotropic, triaxial antenna according to the measurement procedure listed in 32 . We measured total incident E RMS between 0.016 V/m and 0.226 V/m in five rural environments with a linear average of 0.06 V/m and a quadratic average of 0.1 V/m. Joseph et al. 32 measured a median total E RMS value of 0.09 V/m over several rural locations in Belgium, the Netherlands, and Sweden. Bhatt et al. 1 measured an average E RMS value of 0.07 ± 0.04 V/m in rural environments in Belgium. Both previous studies of rural RF-EMF exposure are close to what we found in this manuscript and certainly within the measurement uncertainty of 3 dB on our measurements.
As our RF-EMF exposure measurements near bee hives demonstrate, see Table 2, most of the current RF-EMF exposure is located at frequencies ≤1 GHz. Additionally, Fig. 5 demonstrates that the P abs in all studied Honey bee models is lowest at frequencies ≤1 GHz. This implies that in reality, potential shifts in telecommunication frequencies to higher frequencies might induce even larger increases that the ones estimated in Table 4 since in that analysis an average value over all P abs values ≤3 GHz is assumed.

Strengths and limitations.
This manuscript presents several contributions to the state of the art in the field of RF-EMF exposure assessment of insects. First, to the best of the authors' knowledge, this is the only paper where a numerical RF dosimetry is presented for different developmental stages of honey bees. Second, this is the only study that combined real, in-situ exposure measurements with numerical simulations of RF-EMF exposure of insects in order to estimate a realistic exposure of honey bees. In comparison to our previous study 11 , we considered a broader frequency range from 0.6 GHz up to 120 GHz, which is more in line with the frequencies used in the current telecommunication networks (3 G and 4 G). Finally, this study presents a unique quantification of real-life exposure of honey bees and estimations of how this might change if future frequency shifts in that exposure might occur. A disadvantage of this study is that we did not executed dielectric and thermal measurements in order to obtain dielectric and thermal properties of the studied honey bees. We obtained dielectric properties from literature and were able to execute electromagnetic simulations. We did not perform thermal simulations in this study. Another disadvantage is that we modeled far-field exposure by a limited number of plane waves, while previous studies have shown that a large set of plane waves is necessary to properly model far-field exposure 26 . We did executed a validation of our exposure set up by comparing it with a set of random plane wave exposures and found good correspondence, certainly close to the mean/median. Finally, we used FDTD simulations that are faced with uncertainties 29 and used models that have a limited spatial resolution. This is a disadvantage of any RF-EMF simulation study in comparison to a study that relies on measurements of real insects.

Future research.
Our future research will focus on executing exposure measurements of insects in order to validate the RF-EMF P abs values and the dielectric parameters. Additionally, we would like to execute thermal simulations of honey bees and other insects under RF-EMF exposure. Finally, we aim to work on the development of more insect phantoms, with more spatial accuracy and potentially several independently identified tissues.

conclusions
Exposure of Western Honey Bees (apis mellifera) to radio-frequency (RF) electromagnetic fields was studied using a combination of in-situ exposure measurements near bee hives in Belgium and numerical simulations. The simulations use the finite-difference time-domain technique to determine the electromagnetic fields in and around five honey bee models exposed to plane waves at frequencies from 0.6 GHz up to 120 GHz. These simulations lead to a quantification of the whole-body averaged absorbed radio-frequency power (P abs ) as a function of frequency. The average P abs increases by factors 16 to 121, depending on the considered phantom, when the frequency is increased from 0.6 GHz to 6 GHz for a fixed incident electric field strength. A relatively small decrease in P abs is observed for all studied honey bees between 12 and 120 GHz. RF exposure measurements were executed on ten sites near five different locations with bee hives in Belgium. These measurements resulted in an average total incident RF field strength of 0.06 V/m, which was in excellent agreement with literature. This value was used to assess P abs for those honey bees at those measurement sites. A realistic P abs is estimated to be between 0.1 and 0.7 nW for the studied honey bee models. Assuming that 10% of the incident power density would shift to frequencies higher than 3 GHz would lead to an increase of this absorption between 390-570%. Such a shift in frequencies is expected in future networks.