1 Introduction

Population surveys are important tools for monitoring changes in the size or range of populations over time and thus addressing key questions in ecology (Williams and Brown 2019; Carvalho, de Oliveira et al. 2021), invasion biology (Muirhead, Gray et al. 2008; Phillips, Gillis et al. 2019), and conservation (Didham, Basset et al. 2020; Wagner, Grames et al. 2021). A range of standard methods are available for surveying insect populations, including malaise traps, light traps, pitfall traps, and pan traps (Montgomery, Belitz et al. 2021). Yet, some insects are challenging to survey by these means due to patchy distributions in time or space, low likelihoods of entering traps, or social life histories in which the density of individuals is not readily equated to effective population size (Visscher and Seeley 1982). In such cases, taxon-specific approaches based on knowledge of the target species’ behaviour may provide more efficient, accurate, and cost-effective means of monitoring populations.

Honey bees (Apis spp.) have a key role as pollinators in many natural ecosystems and agro-ecosystems across the world (Tatsuno and Osawa 2016; Hung, Kingston et al. 2018; Baylis, Lichtenberg et al. 2021). The Western honey bee (Apis mellifera) is native to Europe and Africa and has been naturalized across the Americas, Asia, and Australia in recent centuries in association with its use for honey production and crop pollination. The remaining Apis spp. (c. 8 species) are native to Asia, with the tropical forests of South-East Asia exhibiting the highest species diversity (Oldroyd and Wongsiri 2006; Yang et al. 2011). Apis colonies typically comprise a single reproductive female (queen) and thousands of non-reproductive females (workers). During a colony’s reproductive phase, it also produces males (drones) and new queens. Mature drones leave the colony each day at specific times (Koeniger and Wijayagunasekera 1976; Yoshida, Saito et al. 1994) and aggregate in locations known as drone congregation areas (DCAs) to mate with virgin queens on the wing, many metres above the ground (Zmarlicki and Morse 1963). The aggregation behaviour of drones has been harnessed as a means of population survey in Apis mellifera (Utaipanon, Holmes et al. 2019a; Utaipanon, Schaerf et al. 2021). First, several hundred drones are sampled per DCA by luring them into traps held aloft (Williams 1987). Such traps are baited with synthetic Queen Mandibular Pheromone (QMP), a major component of the pheromone bouquet that virgin honey bee queens use to attract drones at DCAs (Butler and Fairey 1964; Boch, Shearer et al. 1975). Second, sampled drones are genotyped and their sibship inferred to estimate the number of local colonies contributing drones to the sampled DCA. Sibship is determined typically by using maximum likelihood methods implemented in programs such as COLONY (Jones and Wang 2010). Finally, the number of contributing colonies is used to infer local colony density, based on an estimate of the typical maximum range over which drones fly from colony to DCA (Utaipanon, Holmes et al. 2019a). This method is advantageous as A. mellifera colonies are often cryptic and difficult to locate, particularly in forests or areas with difficult terrain. Drone sampling, in contrast, allows the detection of many colonies from a local area via a visit to a single site (Utaipanon, Holmes et al. 2019a) and, as drones are actively attracted to lures, may be effective at finding evidence for the presence of colonies, even at low density. However, this approach to population surveys has not yet been applied to other species of Apis.

The Asiatic Hive Bee, A. cerana, is native to Asia, including South-East Asia, China, Japan, and India (Figure 1), and fulfils an important role as a pollinator in native forests across this range (Arbol, 2013). This species is increasingly threatened in much of its native range due to land clearing, climate change, overharvesting, pesticide use (Yasuda, Sakamoto et al. 2017), and competition with non-native managed A. mellifera (Oldroyd and Nanork 2009; Abrol 2013). These combined stressors are proposed to be driving population declines of A. cerana throughout Asia (Abrol 2013). At the same time, A. cerana has recently become invasive in some regions of the Austral-Pacific (Australia, New Guinea and Solomon Islands; Figure 1), where it poses a threat to resident bee species and agriculture via competition and the transfer of pests and pathogens (Anderson, Annand et al. 2010; Koetz 2013b). In both circumstances, there is thus a timely need to investigate new tools for monitoring A. cerana populations.

Figure 1
figure 1

The global distribution of the honey bee Apis cerana. Green indicates the native range of all A. cerana. Red ranges indicate the invasive range of A. cerana in the Austral-Pacific, which includes New Guinea, the Solomon Islands, and Far North Queensland, Australia (adapted from Koetz 2013b).

Full size image

Colony density estimates via DCA sampling rely first and foremost on being able to reliably locate DCAs in the environment. The locations and times of drone flights to DCAs are species-specific (Koeniger and Wijayagunasekera 1976) and sometimes even population-specific (Koeniger and Koeniger 2000). A. cerana has been reported to fly in temperatures above 18 °C, during high humidity, and drones may congregate as low as 2 m above the ground (Punchihewa, Koeniger et al. 1990) or as high as 50 m above the treeline (Fujiwara, Miura et al. 1994), with these differences possibly explained by variation between subspecies (Yang et al. 2011). DCA sampling also relies on a threshold number of drones per colony being recruited to DCAs during the sampling period to be representative of local colony density (Utaipanon, Schaerf et al. 2021). Drones do not fly in all weather conditions, yet the exact conditions that prevent drone flight are poorly understood. Therefore, assessing conditions required for recruitment of drones to DCAs in A. cerana is key to assessing the efficiency and accuracy of drone sampling methods for population surveillance of this species.

In this study, we investigate drone congregation behaviour, and its application to population surveys, in a population of A. cerana in Far North Queensland, Australia. This population arose from an accidental introduction. It was first detected in the port city of Cairns in 2007 (Koetz 2013b) and was likely founded by the arrival by boat of just one or a few swarms from Papua New Guinea, where the species is also invasive (Gloag et al. 2019; Sun et al. 2022). Here, we use a modified version of the William’s Drone Trap, the device commonly used to trap A. mellifera drones, to take advantage of A. cerana’s natural reproductive behaviours for population surveys. In this case, population surveys serve to assess changes in the population’s range and size relevant to its ongoing invasion. Specifically, we (i) locate and describe A. cerana DCAs in this population, (ii) identify the weather conditions in which A. cerana drones congregate, and (iii) identify genotype drones sampled at these DCAs to estimate the density and size of the population during our 5-year study period (2016–2021). We then discuss both the challenges and opportunities for using DCA sampling to monitor populations of A. cerana in both their native and invasive range.

2 Methods

2.1 Drone congregation areas in Australia’s Apis cerana population

Australia’s A. cerana population occurs in a region of Far North Queensland that incorporates rainforest, woodland, farmland, and urban areas. This known occupied region, defined by citizen and researcher reports to Queensland Department of Agriculture and Fisheries (DAF), is estimated at 10,200 km2 in 2022. We searched for and sampled A. cerana DCAs within this invasive range during four field periods: July–August 2016, May–August 2018, July–November 2019, and April 2021. We identified potential A. cerana DCAs for sampling by first inspecting landscapes on Google Earth, where we looked for patches of trees or remnant forest which were accessible by vehicle. We then visually inspected candidate locations on the ground and selected sites where treelines acted as a shelter from wind and thus would allow the use of a balloon. In 2016, we focused on the most accessible sites, which were those in parks, nature reserves, and roadside vegetation in urban areas. In 2018–2021, we revisited successful sites from previous years and also searched locations in natural forest sites near Barron Gorge, Lake Morris Reserve, Wooroonooran National Park, Japoon National Park, Mount Mackay National Park, Djiru National Park, and Daintree National Park. Additionally, in 2018, we attempted to find DCAs and sample drones from the low-density southern range edge, which served to test whether this sampling technique was effective for surveying even where A. cerana abundance was very low. In total, 82 locations were visited over this 4-year period, at which 42 were sampled at least one drone, 30 were sampled at least 10 drones in a single day (these we deemed “confirmed DCAs”), and 11 were sampled ≥ 100 drones across one or more days in the same month (Figure 2A, Supplementary Table 1). In total, 10 DCAs were revisited at least once, including 4 DCAs revisited in all 4 sample years.

Figure 2
figure 2

Sites visited over the course of this study in Apis cerana’s known invasive region in Australia as of 2022 (outline shown in grey; approximately 10,200 km2). A All candidate DCA sites included in this study. Sites where drones were present (≥ 1 drone caught) are shown with larger points and sites where they were absent are shown with smaller points. B DCAs used in this study that were used to calculate population density (all of which had ≥ 100 drones caught; N = 9 DCAs). The colour of the points indicates the number of years these sites were sampled during our study period (2016–2021).

Full size image

We first attempted to catch drones with a traditional William’s Drone Trap, commonly used to catch A. mellifera drones (Figure 3A) (Williams 1987). These traps consist of a net affixed to retractable nylon line which is held aloft in the air using a large helium inflated balloon (90 cm Latex Balloons, Qualatex brand). Attached to the interior and bottom of the net are nylon lines connected to fake queens, which are blackened cigarette filters dosed with 10 µL of E-9-oxo-2-decenoic acid, the main component of queen mandibular pheromone (9-ODA, Intko Supply Ltd.). A. mellifera drones are attracted to the fake queens and fly towards them in an attempt at copulation. At the last second, they realize their mistake and fly upwards towards the hood of the trap from which they are unable to escape. However, after deploying this trap for three consecutive days at each of two separate sites in Cairns (see Supplementary Table 1), it became clear that a modified design would be needed. A. cerana drones were present at these sites and were observed to fly rapidly past the lures at the bottom of the trap. However, they were timid and would not linger near the lures nor did they ever enter the trap. When we reeled in the trap, removed it from the balloon, and instead attached a lure directly to the balloon line, we observed that dozens of drones now approached the lure very closely, even forming “drone comets” around the lure as we moved the balloon. We therefore instead designed a version of the Williams Drone Trap (Figure 3B–C) that allowed us to trap A. cerana drones. This trap was comprised of a retractable nylon line (0.4–1.0 mm) connected to a large (approximately 90 cm) helium-inflated balloon, with 3–5 short (10 cm) heavy gauge (1.0 mm) secondary lines tied to the main balloon line. Each secondary line was tied around a fake queen lure, dosed with the same chemical (9-ODA) as for A. mellifera drones. A 9-ODA is similarly the major component of A. cerana queen mandibular pheromone (Plettner, Otis et al. 1997), and A. cerana drones are strongly attracted to it (Koeniger and Koeniger 2000). Each of the smaller lines was distributed approximately 20 cm from each other, the topmost being approximately 1 m below the balloon. We then coated the lures, tethers, and up to 1m of the nylon line below the last lure of the trap with flypaper glue (The Buzz: Fly Papers Insect Trap—73,018). We then raised the balloon to the height of the treeline, such that the queen lures lay slightly below the treeline. Once our modified Williams Balloon Trap had been raised at a candidate site, we moved along the treelines until we observed A. cerana drones. We then continued to move along the treeline and vary the balloon height until we located the position at which the greatest number of drones showed interest in the lures (this occurred between 2 m and the height of the treeline). Unlike A. mellifera, we observed that A. cerana drones do not remain airborne during mating flights but instead dart out from nearby vegetation. During this time, we assume that drones rest on leaves and make periodic flights to inspect lures, although we could not observe this behaviour at that height. When A. cerana drones approached the queen lures in efforts to copulate, they touched the line and became stuck on the glue that coated the trap (Supp Video 1 & Supp Video 2). We then reeled the line in and preserved drones in 100% ethanol. Additionally, in 2021, for some DCAs, we measured the approximate height of the tallest tree in the treeline that surrounds a DCA, and the tallest tree within eyesight (~ 100 m) of the DCA. We measured tree height by allowing the bottom of the balloon to reach the approximate tree height, reeling it down and then measuring the length of the line.

Figure 3
figure 3

Different drone traps used to catch A. mellifera (A; Williams 1987) and A. cerana (B; this study) respectively. A A traditional Williams Drone trap, wherein drones are attracted to fake queen lures and find themselves caught within the net as they try to escape. B Our modified Williams Drone trap, where drones are attracted to fake queen lures and find themselves caught on the line which is coated in flypaper glue. C An example of balloon height from a DCA used in this study. Sample videos on catching drones using this method can be found here (Supp Video 1) and here (Supp Video 2).

Full size image

Trapping took place between 1 and 4 pm Australian Eastern Daylight Time (AEDT). We aimed to sample at least 100 drones per DCA per afternoon. If a DCA had a high density of drones, it was possible to secure 15–20 drones on the line each time it was reeled in. Where drone density was lower, we sampled 1–5 at a time. Glue was reapplied every 10 min during active trapping. Where possible, sites with fewer than 100 drones caught in a single day were repeatedly sampled over multiple days to increase sample size. As bees do not fly in the rain (Reyes, Crauser et al. 2019), we did not attempt to sample on any days with precipitation. We recorded the size (m2) of high abundance DCAs (≥ 100 drones caught in one field trip/year; N = 11) and describe these in the results. Lastly, we describe the physical factors (size, surrounding treelines and drone flight pathways) of three DCAs with the highest sample sizes over our study period.

2.2 Abiotic factors that influence drone recruitment to DCAs

During our efforts to locate and sample A. cerana DCAs, we noticed that the number of drones present at a site varied between days, even when skies were clear and without rain. To better understand how weather affects drone recruitment to DCAs, we collected climate data for all days we attempted to trap drones at a subset of our DCAs that had high drone abundance (≥ 100 drones caught in one field trip/year; N = 11 DCAs; Supplementary Table 1). For each DCA, we determined the closest Australian Government Bureau of Meteorology (BOM) weather station that recorded weather conditions (wind speed (km/h), temperature (C°), and relative humidity (%)) by the minute. Distances between DCAs and BOM weather stations ranged from 1 to 16 km (mean = 7 km). For each DCA and sampling day, we recorded drone presence (at least one drone trapped), drone absence (no drones trapped), and (if present) number of trapped drones and also found the average wind speed, temperature, and relative humidity between 12:00 and 16:00 AEDT. We then performed (i) logistic regressions on drone presence/absence for all weather conditions using the glmer function in R (4.0.4), using site as a random variable, and (ii) for variables that returned a significant logistic regression, a linear regression of the variable and drone number using the lmer function in R (4.0.4), using site as a random variable.

In addition, we collected data on local weather conditions at the DCAs at the time of sampling for DCAs visited in 2021 only (8 DCAs with some revisited, N = 14; Supplementary Table 1). We did this to better understand whether drones selected aggregation sites that differed in abiotic properties to the regional average (e.g. more sheltered from wind, or warmer, than average values measured at weather stations). We collected local weather data in 2021 by sampling wind speed (km/h), temperature (C°), and relative humidity (%) over a 1-min interval every 15 min from our arrival (earliest 13:00 and latest 14:00) until 15 min after the last drone was seen (earliest 14:45 and latest 15:15) using an Ozito digital handheld humidity and temperature meter (ODTH-100) and Meteos/Atmos Anemometer (IC-AME-00). We then assessed whether local and regional weather conditions were correlated with one another by performing a linear regression using the glm function in R (4.0.4). If conditions were not correlated with one another for a given variable, we then performed a paired t-test using the t.test function in R (4.0.4) to understand how local weather conditions within DCAs varied relative to the regional weather.

2.3 Estimating population size and density

We used the drones we sampled from DCAs to estimate the density of colonies in Australia’s invasive population of A. cerana, following the protocol used previously for A. mellifera (Utaipanon, Holmes et al. 2019a). For this analysis, we considered samples from 11 DCAs where we had sampled at least 100 drones. This included 4 DCAs that were represented over multiple sampling years (N = 2 for 2 years; N = 3 for 3 years; N = 19 total DCAs used for this analysis; Supplementary Table 1). We genotyped all drones caught at these DCAs at eight microsatellite loci (Ac27, Ac32, Ac35, A107, Ac1, Ac26, B124, and Ac3 (Takahashi, Shimizu et al. 2009)). As Australia’s population of A. cerana experienced a genetic bottleneck during founding, genetic diversity at these microsatellite loci is low compared to native range populations (Gloag, Ding et al. 2017). Loci had an average of 3.5 alleles per locus (range: 2–5), and missing call rates averaged 3.9% (range: 0.4–19%, although the second highest rate was 4.4%; Supplementary Table 2–3). We extracted DNA from one hind leg per drone using a Chelex protocol (Walsh et al. 1991) and amplified DNA in 5 μl reactions (1 × reaction buffer, 2.5 mM MgCl2, 0.16 mM dNTP mixture, 0.32–0.8 μM of fluorescent dye-labelled primers, Sigma-Aldrich, St Louis, MO, USA, 0.4 units Taq polymerase and 1 μl extracted DNA) and visualized products using an ABI3130xl genetic analyzer (Applied Biosystems, Foster City, CA, USA). For each DCA collection, we then performed a pedigree reconstruction using drone genotypes in the program COLONY (Jones and Wang 2010). In honey bees, provided that a queen is present, almost all drones are sons of the queen (Oldroyd, Halling et al. 2001). Thus, every mother that was inferred from the pedigree analysis indicates a colony in the vicinity of the DCA; that is, the number of mothers is equal to the number of colonies (Nc) represented at a DCA. However, A. cerana workers are known to occasionally activate their ovaries and produce haploid drone offspring, the presence of which at DCAs interfere with the pedigree reconstruction of the COLONY program. Such worker-borne drones are typically laid in worker cells and are smaller on average than queenright drones. Therefore, we measured head sizes for all drones and excluded drones with heads smaller than 3.2 mm from the pedigree analysis (as drones this size are highly likely to be worker-laid (Gloag, Christie et al. 2019; Utaipanon, Holmes et al. 2019b)). We also excluded DCAs from further analysis if the removal of small drones resulted in a drone count of < 100 (this removed three sites, resulting in N = 16 DCAs used for the final analysis).

Unfortunately, pedigree analysis of haploid (male) honey bees can be inaccurate at low sample sizes (Utaipanon, Schaerf et al. 2021; Williamson, Groom et al. 2022). If a colony is represented by single drone at a DCA, there is a high probability that it will be clustered into a family with another singleton, rather than each singleton being correctly assigned to their own unique mother. This is because the algorithm implemented in COLONY searches for the minimum number of families that can be explained by the genotypes analysed. Also, when a colony is represented by just a few drones, the algorithm may split them into two sibships (Williamson, Groom et al. 2022) or re-assign them into genotypically similar families. This error is known as the misassignment (Williamson, Groom et al. 2022) or non-detection error (Utaipanon, Schaerf et al. 2021), as drones present in the sample are misassigned and families go undetected. Previous research in A. mellifera indicates that the misassignment error only becomes significant when the average number of males (DC) per sibship (or the number of local colonies, NC) is < 6, and excluding sampled sites that fall below this threshold (i.e. Dc/Nc = 6) results in a more robust estimate of the real number of colonies surrounding a DCA (Utaipanon, Schaerf et al. 2021). Therefore, we only included a DCA in our final density estimates if the ratio of drones to colonies (Dc/Nc) for the DCA was above 6 (N = 14 from 9 DCAs over several years, Supplementary Table 4).

To estimate colony density surrounding each DCA, we divided Nc by an estimate of the catchment area that surrounds a DCA. We consider this area to be a circle surrounding the DCA equal to the maximum distance drones may fly to a DCA from their colony of origin. This distance is not known for A. cerana, but in A. mellifera, it has been estimated to be 3.75 km (catchment = 44.2 km2 (Utaipanon, Holmes et al. 2019a)). However, flight capacity in bees is related to body size (Greenleaf, Williams et al. 2007) and A. cerana drones are smaller than A. mellifera drones (intertegular distance: 3.2 ± 0.1 vs 4.3 ± 0.1 mm respectively (Streinzer, Brockmann et al. 2013)). Using the equation of Greenleaf et al. (2007), we can approximate a flight range for A. cerana drones, where range = YOMb, and YO is a constant, M is intertegular distance, and b = 2.7. Specifically, we considered that the flight range of A. cerana drones is determined by body size proportionally to the flight range of their closest relative, A. mellifera (RangeCerana = RangeMellifera × YOMCeranab/YOMMelliferab). Using this relationship between flight range and body size, A. cerana drone flight distance should be approximately 1.70 km (catchment = 9.0 km2). Given the uncertainty regarding true drone flight distance of A. cerana, and the importance of this variable for colony density estimates, we therefore consider two possible values for A. cerana drone flight distance: 3.75 km (i.e. the flight distance for the larger-bodied A. mellifera and thus a conservative estimate for A. cerana) and 1.70 km (the estimate obtained from a power relationship with body size according to Greenleaf et al. 2007). To estimate the number of colonies in this population (population size), we multiply our density estimates by the current area the population spans (all localities that A. cerana has been confirmed as invasive as of 2022), which is approximately 10,200 km2.

Finally, as our methods rely on colonies presenting drones at local DCAs, we will underestimate population size and density when there are local colonies that are not currently producing drones, or producing them in very small numbers. This error rate is sometimes called the non-sampling rate (Utaipanon, Schaerf et al. 2021; Williamson, Groom et al. 2022), as colonies that do not produce drones cannot be sampled. Some studies that investigate A. mellifera regard this error rate as unimportant (Utaipanon, Schaerf et al. 2021), on the basis that all strong A. mellifera colonies produce drones seasonally at the same times of year and colonies that are not strong enough to do so are unlikely to be relevant to agricultural processes such as pollination. However, this error clearly cannot be ignored for A. cerana in Australia and other tropical Apis species where reproductive cycles are highly variable across a year, and many colonies will not be in the drone production phase at any given time (Oldroyd and Wongsiri 2006; Koetz 2013a). Therefore, we multiply our estimates by a factor of 1/α, where α is the percentage of colonies that produce drones at any given time, to account for non-sampling error. In this population, we estimate that 33% of colonies are producing drones at any given point in time (based on sampling of 15 drone-producing nests from 45 total nests collected July–September 2015–2016; Gloag et al 2019). However, as this estimate of non-sampling rate is based on low sample sizes and we could not explicitly test this rate during all years that we sampled drones, we chose a non-sampling rate of 50% (α = 0.5) to provide a minimum estimate of population size and colony density, and a non-sampling rate of 33% (α = 0.33) for our upper estimate.

3 Results

3.1 Drone congregation areas in Australia’s Apis cerana population

We identified 30 DCAs in this population of A. cerana (defined as sites where we trapped at least 10 drones in 1 day). At the three DCAs in this study where we observed the highest drone activity, we sampled ≥ 250 drones in a single day (Figure 4: A, N = 271; B, N = 447; C, N = 414). We revisited 6 sites over multiple years, all of which remained active DCAs. This included four sites that were active across all 5 years of our study (2016–2020). We found that DCAs were small areas of open sky edged by trees on two or three sides (Figure 4). The average area of a DCA was 430 m2, but we found this varied between 100 and 1800 m2 and that 70% of sites had areas less than 300 m2. The trees surrounding a DCA are usually consistent in height and regular in placement. Canopy heights of the trees surrounding DCAs have been previously determined to be 4–8 m (Punchihewa, Koeniger et al. 1990), but within this study (for DCAs sampled in 2021), we found that the highest canopy height of a DCA varied between 11 and 24 m and the height of the highest tree within the local area varied between 18 and 28 m. Notably, while all DCAs that we detected shared these general characteristics, many locations that we sampled that shared these characteristics did not harbour drones, despite these sites being in areas where colony density was known to be high.

Figure 4
figure 4

Three example Apis cerana DCAs in Australia. A–C Approximate height at which drones were caught at these DCAs (shaded regions). D–F Aerial views of the same DCAs, labelled with approximate area (m2; dotted lines) and possible routes along treelines that drones might follow to reach these DCAs (white lines).

Full size image

3.2 Abiotic factors that influence drone recruitment to DCAs

Based on data from the nearest weather stations, the only significant predictor of drone presence at DCAs was temperature (β = 0.60 ± 0.17, z = 3.44, p = 0.0006), where higher temperatures resulted in a higher chance of drone presence. Notably, all peak trapping days (where we trapped ≥ 100 drones in a single day) occurred when temperatures were greater than 24 °C, and days below this threshold had a lower number of drones caught (mean = 7 drones, range: 0–45, Figure 5). We found no significant relationship, however, between drone presence at a DCA and wind speed (β = – 0.07 ± 0.04, z = – 1.66, p = 0.098) or humidity (β = 0.00 ± 0.03, z = 1.33, p = 0.89). And although temperature had a significant impact on drone presence at DCAs, we found no such relationship between temperature and the actual number of drones we caught at these sites (β = 6.87 ± 3.67, t = 1.87, p = 0.064).

Figure 5
figure 5

The relationship between abiotic factors (as measured by nearest weather stations) and drone recruitment to DCAs in Australia’s Apis cerana population. A Average temperature of days on which we trapped < 100 drones vs. ≥ 100 drones, where the latter was considered a particularly good trapping day. B Relationship between average temperature and drone recruitment, with dotted line indicating the lowest temperature at which at least 100 drones were sampled at any site (24 °C). C Relationship between average wind speed and drone recruitment, with dotted line indicating the highest recorded local wind speed of interest at a DCA site (8 km/hr), which was below all regional wind speeds recorded at the nearest weather stations. D Relationship between average humidity (%) and drone recruitment. Blue lines in all cases indicate linear regressions with 95% C.I.s (grey).

Full size image

When we considered the relationship between weather data from stations and local measurements taken at the DCA site, we found that while regional and local temperature were correlated (β = 0.94 ± 0.31, t = 3.01, p = 0.011), our local measurements were consistently 2 °C warmer (Supplementary Table 5). Neither local wind speed (β = 0.01 ± 0.12, t = 0.12, p = 0.91) nor local humidity (β = 0.53 ± 0.25, t = 2.16, p = 0.052) were correlated with conditions measured from weather stations. Moreover, local conditions at DCAs were significantly less windy (on average 14.6 km/h less; t[13] = – 10.8, p < 0.001) and more humid (on average 5.5% more; t[13] = 3.41, p = 0.005) than those measured by weather stations (Supplementary Table 5). Indeed, the local wind speeds we measured at DCAs were always < 8 km/h, while 8 km/h was the minimum of the average wind speed measured by weather stations on our collection days (this ranged from 8 to 36 km/h, mean = 20 km/h); that is, local wind speeds at DCAs were always far lower than regional wind speeds on the same sample day.

3.3 Estimating population size and density

Based on sibship analysis in COLONY, an average of 24 ± 6 colonies contributed drones per DCA (n = 14, range: 17 37). While different years are represented unevenly (2016 = 8, 2018 = 4, 2019/2021 = 2), we find no significant difference between years in the colony numbers represented per DCA (F2,11 = 0.364, p = 0.703).

Using an estimate of drone flight distance of 3.75 km (obtained from A. mellifera) and a colony non-sampling rate of 50% (i.e. proportion of all colonies in the area producing drones), we reach a conservative colony density estimate of 1.1 ± 0.3 colonies/km2 for this population. Assuming an even distribution of colonies across the known range of A. cerana in Australia, this results in a minimum estimated population size of 11,000 ± 3000 colonies. If we instead estimate drone flight distance to be 1.70 km in A. cerana (based on its smaller body size than A. mellifera) and assume a non-sampling rate of 33% (best estimate based on available data from this population; Gloag et al. 2019), then colony density is estimated to be 8.1 ± 2.0 colonies/km2, and population size is estimated at 83,000 ± 21,000 colonies.

4 Discussion

We located drone congregation areas in Australia’s invasive population of A. cerana, sampled drones from these sites using QMP-lures attached to helium balloons, and inferred colony density estimates and population size by genotyping drones and inferring the number of colonies contributing to each DCA. We conclude that a drone-sampling approach to population surveys in A. cerana is feasible and that the method presents both advantages and challenges.

A key requirement of honey bee population surveys from DCAs is that the DCAs themselves can be located. We found that we were able to locate drone congregation areas in Australia’s A. cerana population fairly readily by searching landscapes for sites with small, sheltered areas of open sky adjacent to tree lines. This finding is consistent with studies from native range populations, wherein A. cerana drones prefer DCAs that are more enclosed (Punchihewa, Koeniger et al. 1990) than the large open areas preferred by A. mellifera (Loper and Taylor Jr 1987). In our study, the typical A. cerana DCA was < 300 m2, whereas A. mellifera DCAs average 1600 m2 (Loper and Taylor Jr 1987). Drone congregation times and locations between co-occurring Apis species are expected to diverge in response to selection against interspecific mating (Koeniger and Koeniger 1991, 2000). For example, A. cerana males in Sri Lanka fly at different times, and in different locations, to those of sympatric A. florea and A. dorsata (Koeniger and Wijayagunasekera 1976). However, given that A. cerana and A. mellifera have been brought into sympatry in Australia only recently via human-induced range shifts (Remnant, Koetz et al. 2014), different DCA locations in this case may principally reflect differences in the native environments of the two species (i.e. rainforest versus open woodland) rather than selection against interspecific mating. Variation in local environment may also explain variation in the DCA characteristics described for different A. cerana subspecies. Drones of A. cerana japonica (locality: Japan) have been reported to prefer congregating 30–70 m above the ground (Fujiwara, Miura et al. 1994), while A. cerana (locality: Sri Lanka) are reported to prefer 2–8 m above the ground (Punchihewa, Koeniger et al. 1990). In the Australian invasive population (derived from the native populations of Java, Indonesia; Sun, Chapman et al. 2022), we observed that drones flew close to the top of the treeline, anywhere from 2 to 24 m above the ground.

Notably, some sites that we visited met our physical criteria for A. cerana DCAs and were in areas where A. cerana density was known to be high, yet we failed to sample drones at the site. Therefore, either there are some required features for A. cerana DCAs that remain unknown to us or when multiple suitable DCA sites occur in a local area, the one at which males actually congregate is stochastically decided via drone and/or queen critical mass (e.g. some drones may choose one site by chance and thereafter attract all further drones to that site). Additionally, we cannot rule out that A. cerana drones in Australia also sometimes congregate at different types of sites, given that we principally searched only locations that were both sheltered and accessible enough for us to deploy the balloon (i.e. we did not search systematically every possible kind of site). Nevertheless, because some of the sites we did find regularly attracted upwards of 100–400 drones and were stable across multiple years, we can conclude that these types of DCAs are suitable for population surveys. Although many of our sampled DCAs used for population density were in disturbed habitats near human habitation (where access was easiest), we were able to also sample drones from four sites in the forested areas of National Parks using the same search strategy, albeit at lower sample sizes.

Surveys of Apis spp., particularly A. mellifera (Loper, Wolf et al. 1992), based on DCA sampling are sensitive to the weather. We found that warm temperatures are a prerequisite for the recruitment of drones to DCAs in Australia’s invasive A. cerana population, with zero or few drones captured at temperatures below 24 °C. Drone production and swarming occur year round in tropical A. cerana populations, and this minimum flight temperature is regularly reached even during the winter months in our study area (June–August). Yet, it seems likely that the minimum temperature required for mating flights may vary significantly across the native range of A. cerana given its large area, spanning both tropical Asia and temperate regions of Japan, India, Pakistan, and China (Hepburn and Radloff 2011). Furthermore, variables other than temperature must be involved with day to day drone recruitment in Australia, as we sometimes found that drones were not recruited to known DCAs even at high temperatures (≥ 24 °C). Wind speeds in particular are likely to be important, given that we found local wind speeds at DCAs were always lower than wind speeds recorded from nearby weather stations, consistent with the idea that drones select aggregation sites protected from excessive wind. In this way, A. cerana drone behaviour would also be consistent with that of A. mellifera, for which high winds result in fewer drones leaving the colony (Reyes, Crauser et al. 2019).

A. cerana population surveys via DCA sampling can have two broad goals. The first is the same goal previously proposed for A. mellifera: estimating local colony density and population size. While our study shows that this is achievable in A. cerana, several aspects of Asian honey bee biology are currently poorly known and present noteworthy challenges to this goal. These include the following: (1) the typical drone flight distance of males; (2) the non-detection error introduced through pedigree analysis; (3) the non-sampling error (i.e. the typical proportion of colonies in the area that are not producing drones) and; (4) the heterogeneity of colony density across space. In this study, we take various estimates and approaches for these four measures to arrive at lower and upper estimates of colony density. First, we use two estimates of A. cerana drone flight distance: one which assumes a flight distance similar to the larger-bodied A. mellifera (3.75 km, the only Apis species for which such estimates have yet been made; Utaipanon, Schaerf et al. 2019) and one which assumes a drone flight distance smaller than that of A. mellifera, approximated from body size (1.70 km; Greenleaf et al 2007). The assumed value for typical drone flight distance has a large impact on the estimate for colony density, and data on the true flight range of A. cerana drones will be essential for making more refined estimates of colony density in the future. Second, we only include DCAs which have a minimum non-detection error based on the minimum required ratio of drones to mothers (DC/NC) obtained from studies in A. mellifera (Utaipanon, Schaerf et al. 2021). Unfortunately, Australia’s invasive A. cerana population presents a further challenge for sibship analysis due to low genetic diversity (Gloag et al. 2017) and the presence of drones reared in queenless colonies that are the sons of workers (arising from high rates of queenlessness in this population; Gloag, Christie et al. 2019)). These factors could bias the number of mothers calculated by COLONY for each DCA, resulting in lower or higher colony density respectively. In this study, we chose to exclude the smallest drones on the basis that they may have been worker borne, but doing so may have biased our sibship analysis towards a lower estimate of colony density. Moreover, the low polymorphism at the markers we used may also lead us to underestimate colony numbers via sibship misassignment errors in COLONY. This latter issue, however, is thankfully rather unique to Australia’s invasive population (due to its very severe founding bottleneck) and should not present a problem for researchers employing this technique in native range populations. Third, we assumed colony non-sampling rates (α) of either 50% (a conservative estimate) or 33% (the rate indicated for Australia’s population based on available data; Gloag et al. 2019). Whether this range of α is suitable for estimates in all populations of tropical A. cerana however is unclear. Ideally, some data on the mean and variance of this parameter would be available for a given population before attempting to estimate colony density via the drone sampling approach. Finally, in the absence of data from every corner of the population (Figure 2), we assume a homogenous colony density across the species’ invasive range (currently just 10,200 km2). This final assumption would artificially inflate our estimates of total population size, if non-sampled areas in fact tended to have lower population densities than those that we sampled. However, given that we found population densities were consistent over large distances (Figure 2) and did not vary between the urban and non-urban sites in our dataset, we consider that high spatial heterogeneity of colony density is unlikely (with the exception of the extreme range edges, where density is lower). When applying this survey method in other A. cerana populations, knowledge of the local A. cerana will help guide the best approach to estimating drone flight distance, non-sampling error and spatial heterogeneity, as well as the choice of number and type of molecular markers for sibship analysis.

A secondary goal of DCA sampling is to assess the presence of Apis spp. in areas where density is expected to be very low. Such low-density areas may include recently-established invasive populations, or the range edges of expanding or contracting populations. Here, DCA sampling can help to detect the presence or absence of Apis. spp., including A. cerana, in locations where direct detection of colonies may be extremely difficult due to their low density and cryptic nature. In fact, the use of this technique has already proved useful in this way for invasive A. cerana in Australia: The modified William's Drone Traps was used during the biosecurity response to a secondary incursion of A. cerana in Townsville, Queensland in 2016, when searching for colonies to eradicate (ultimately this population failed to establish (DAF 2019)).

In conclusion, we propose that modified William’s Drone Traps used at A. cerana DCAs can be part of a toolkit of methods used to survey and manage A. cerana in both its native and invasive ranges, as it is for A. mellifera (Beaurepaire, Kraus et al. 2014). Through this method, we estimate the size of Australia’s A. cerana population in the range of 11,000–83,000 colonies, a density of around 1–8 colonies per square kilometre. As research advances on the biology of A. cerana drones and drone production, and the drones of other Apis species native to Asia, the accuracy of colony density and regional population size estimates via this method will increase.