Global invasion risk of Apocephalus borealis, a honey bee parasitoid

Apocephalus borealis is a parasitoid of hymenopterans native to North America that also attacks introduced honey bees (Apis mellifera). Parasitism by this species has been associated with infested bees absconding the hive and dying outside. The flies can also harbour viral infections and nosematosis. Recently, nucleotide sequences identical to A. borealis were reported from bulk screenings of honey bees from Belgium and South Korea, although no adult flies have been collected. To predict the potential invasion risk of A. borealis across the world, we constructed a MaxEnt species distribution model based on occurrence data from North America submitted to the citizen science project ZomBee Watch (zombeewatch.org) and from museum specimens. The results have shown that extensive parts of Europe, the Mediterranean Basin, Asia Minor, southern Africa, eastern Asia, Australasia, and North and South America have high degrees of climatic suitability for invasion, suggesting that the fly could establish in these regions. The potential invasion range is expected to stay similar under different climate change scenarios. We discuss the status of A. borealis as an invasive species and measures that can be taken to reduce the risk of its introduction outside of North America. Our results highlight A. borealis as a potential threat to honey bee health worldwide that requires urgent attention of international veterinary bodies to prevent its spread.


INTRODUCTION
Apocephalus borealis, or the zombie fly as it became colloquially known, is a phorid fly (Diptera: Phoridae) native to North America that parasitizes diverse hymenopterans, namely, the yellowjacket wasps Vespula arenaria, Vespula pensylvanica, and Vespula vulgaris (Ennik 1973), the carpenter bee Xylocopa varipuncta (Brown 2017), and several species of bumble bees (Brown 1993;Disney 1994). In late 2009, the fly was first recorded to infest adults of the introduced honey bees (Apis mellifera), documenting a recent expansion of its host range (Core et al. 2012). Females (Fig. 1) oviposit into adult bees and the larvae feed on the thoracic muscles and possibly parts of the nervous system of their still living host. Parasitized adult bees appear disoriented, abscond their nest at night, and die outside (Core et al. 2012). These symptoms resemble those of colony collapse disorder (CCD), a major source of bee colony mortality in North America over the last two decades (vanEngelsdorp et al. 2009;Ellis et al. 2010;Lee et al. 2015). Larvae emerge from between the head and thorax of their dead host to pupariate, up to 13 larvae can emerge from a single host in field conditions. Core et al. (2012) demonstrated that some phorid larvae and adults harbour Nosema ceranae and the deformed wing virus (DWV), suggesting that zombie flies may transmit bee pathogens linked to colony mortality (de Miranda and Genersch 2010;Fries 2010;Steinhauer et al. 2018). In bumble bees, A. borealis infestation shortens adult lifespan by up to 70% (Otterstatter et al. 2002), but it is currently difficult to draw firm conclusions about the effects of phorid parasitism on bumblebee and honey bee populations (Cohen et al. 2017).
The fly is a relatively abundant parasitoid in some parts of North America. For example, in sampled Californian apiaries up to 17% of tested foragers were infested (Core et al. 2012). At present, virtually nothing is known about the ecological requirements of A. borealis, fly survival rates on different hosts, its dispersal ability, or the extent to which it may contribute to colony losses. Besides, these worrying gaps in our knowledge of this overlooked but potentially devastating bee parasitoid, it is even more alarming that in recent years its potential occurrence has been reported from outside of its native range. In 2013, a molecular study of honey bee workers from over 360 colonies stationed across Belgium revealed the presence of amplicons showing a 100% similarity to A. borealis in 31% of the studied samples (Ravoet et al. 2013  To predict the geographic range of the parasitoid, understand its ecological requirements, and identify high-risk areas prone to invasion, we built a maximum entropy presence-only species distribution model (SDM) using six bioclimatic variables. Occurrence data were obtained from publications, museum specimens, and ZomBee Watch, a citizen science project encouraging members of the public and beekeepers to monitor the presence of A. borealis in honey bees (zombiewatch.org; Hafernik 2020). We furthermore tested the impacts of different climate change scenarios on the potential invasion range of A. borealis.

Apocephalus borealis occurrence data
Presence data for A. borealis in North America were obtained from ZomBee Watch (zombiewatch.org), a citizen science project that aims to monitor the distribution of the phorid, in January 2021. The ZomBee Watch World Wide Web site was designed and built jointly by members of the San Francisco State University Departments of Biology and Computer Science using modern Agile and user-centred Design methodology. The project, open to members of the public and beekeepers, provides information on how to recognise parasitized bees and construct light traps for monitoring. Participants are asked to upload photos of bees, hold them in containers, and then upload photos of any eventual fly pupae and adult flies that emerge. In addition to geolocation and occurrence information, when needed to confirm identification participants are asked to submit preserved samples of flies or pupae for identification. Further occurrence records were obtained from the online Phorid Catalog (PCAT; https:// www. phorid. net/ pcat/ index. php), based on museum specimens. Additional occurrence data from North America after 1990 were obtained from entomological papers (Otterstatter et al. 2002;Brown 2017). In total, 123 occurrence points were amassed (Fig. 2).

Climatic variables
Grid-referenced continuous bioclimatic variables were downloaded from CliMond version 1.2 (Kriticos et al. 2012), a global database of 19 core monthly bioclimatic data layers averaged between 1960 and 1990, at 10' spatial resolution. These layers were converted to ASCII format using QGIS 3.19 (QGIS Development Team 2020). After initial experimentation and assessment of collinearity, six biologically realistic environmental variables were used as predictors in our model: isothermality (BIO3), temperature seasonality (BIO4), maximum temperature of the warmest month (BIO5), annual precipitation (BIO12), wettest month precipitation (BIO13), and driest month precipitation (BIO14). These variables were chosen based on the fly's known biology and ecological requirements, largely mirroring those used in previous studies modelling the distribution of dipterans (Petersen 2013; Carvalho et al. 2015;Ibáñez-Justicia et al. 2020;Wan et al. 2020). To predict the impacts of climate change on the biogeographical range availability of the zombie fly, we used predicted climate values for 2030 and 2050 derived from two global climate models (GCMs): CSIRO-MK3.0, and MIROC-H. For each model and timeframe, two IPCC IV Special Report on Emissions Scenarios (SRESs) were explored: A1B (medium emissions) and A2 (high emissions) (Nakicenovic and Swart 2000).

Habitat suitability modelling
Habitat suitability for A. borealis for present day and future predicted environmental conditions was computed using the presence-only maximum entropy modelling software (MaxEnt) version 3.4.1 (Phillips et al. 2006(Phillips et al. , 2017. Max-Ent is an SDM algorithm that employs maximum entropy to predict habitat suitability by matching the distribution of occurrences to environmental data and comparing them to background points to estimate the thresholds of species' tolerances and extrapolate species ranges beyond their known distribution (Elith and Leathwick 2009). Models in MaxEnt were built using 123 presence-only occurrence points. Since correcting for geographical sampling bias has been shown to improve the predictive performance of MaxEnt (Syfert et al. 2013), we utilised the kernel density estimator (kde2d) function of the MASS package in R version 4.0.2 (R Core Team 2020) to prepare a raster bias file for inclusion in the analysis (Venables and Ripley 2002). To select MaxEnt settings that balance model fit and predictive ability we used the ENMevaluate function in the ENMeval package in R software, which executes MaxEnt across a range of settings and returns a data file comparing model performance (Muscarella et al. 2014). Based on the ENMevaluate output, MaxEnt was run with linear, quadratic, and hinge features, and with a regularisation multiplier of 3.0. The default "cross-validation" setting with 10,000 randomly selected background points across the study area was used. The final model was computed from an average of 10 runs.

Model performance evaluation
For model evaluation, the area under the receiver operator curve (AUC) parameter of the receiver operator characteristic (ROC) was computed. The AUC can be interpreted as a measure of the probability of correctly predicting species presence in a given geographical cell. An AUC of 0.5 indicates that the model is no better than chance in predicting species distribution, while an AUC of 1.0 represents a perfect fit with no false negatives (Fielding and Bell, 1997;Flø et al. 2015). In addition, analyses of variable contributions and jackknife analyses (Tukey 1958) were performed to assess the contribution of each individual bioclimatic variable to the model output.

Niche overlap testing
Distribution records for A. borealis in Northern America are notably disjunct, divided roughly into a western (99 occurrences) and an eastern cluster (24 occurrences) at the 100th meridian. To test if the eastern and western populations occupy different climatic niches, we measured niche overlap using Schoener's D, which visualises niches in 2D climatic space using a principal component analysis (Schoener 1970), correcting for sampling bias using a smooth kernel density function as described by Broennimann et al. (2012). Schoener's D index ranges from 0 to 1 to describe no to complete niche overlap.
Two measures of niche overlap were employed. The niche equivalency test examines if two niches are indistinguishable from each other (Graham et al. 2004;Warren et al. 2008). Occurrence points from both ranges are pooled together and split into two sets, while maintaining the actual number of occurrences in each range, and the niche overlap statistic D is computed. The process is repeated 100 times. If the value of observed niche overlap, D, falls within the density of 95% of the simulated D values, the null hypothesis of niche equivalency cannot be rejected (Broennimann et al. 2012). On the other hand, the niche similarity test compares if the two ranges are more similar to each other than expected by chance (Broennimann et al. 2012). The density of occurrences in one range is randomly shifted and the overlap of the simulated and observed niches is calculated. As before, the procedure is repeated 100 times.
A p-value > 0.05 indicates that niches are less similar than expected by chance. All analyses were conducted in the ecospat package (Cola et al. 2017) in R. The code as well as the resultant niche equivalency and similarity plots are available from Mendeley Data (https:// doi. org/ 10. 17632/ wvhcm wdwxk.2).

Current distribution range of the zombie fly
Based on occurrence records compiled from ZomBee Watch, the Phorid Catalog, and published references, A. borealis occurs in Canada and the USA. The species extends to Alaska, and the fly has also been recorded from Alberta, British Columbia, New Brunswick, Nova Scotia, Ontario, and Québec in Canada. The records of specimens are concentrated in two main groups, with numerous records from the East and West Coasts but only a few isolated occurrences in mid-North America. This pattern is possibly an artifact of collecting effort by entomologists, as museum collections were not exhaustively queried for further records. Regardless, the niche equivalence test yielded a non-significant result (p = 0.36), failing to reject the null hypothesis. Moreover, the niche similarity test showed that the populations were more similar to each other than expected by chance (p = 0.25). This indicates that the eastern North American zombie fly populations do not occupy a climatically distinct niche from their counterparts in the West Coast.

Predicted distribution of Apocephalus borealis under current climate conditions
MaxEnt distribution modelling showed a high predicted habitat suitability of A. borealis across most of western, southern, and central Europe, parts of the Mediterranean Basin, Asia Minor, the Caucasus, extensive areas of central and southern Africa, eastern Asia, coastal areas of Australia, New Zealand, and coastal and mountainous regions of North and South America (Fig. 3). In addition, habitat suitability was high on islands of the North Atlantic Ocean, Madagascar, New Caledonia, and the Aleutian Islands.
Within North America, the prediction of the SDM obviously corresponds well to the known distribution of the fly, with areas of high habitat suitability found in coastal regions. This is unlikely to represent a modelling artifact since at least in some parts of its distribution range A. borealis seemingly prefers coastal areas. In southern California, the zombie fly is only found within about one kilometre from the ocean, while hundreds of Malaise trap located further inland yielded none (B. Brown, personal observation). However, this may not apply throughout the entire distribution range of A. borealis, where collecting efforts have been mainly concentrated around the coast. For example, in the San Francisco Bay Area, ZomBee Watch has recorded positive samples in arid regions several kilometres East of the Bay.
Another caveat is the likely insufficient known host range of the fly. As the presence of suitable hosts is an important predictor of habitat suitability, future studies of host specialization are needed in order to fine-tune ecological niche models. Understanding the degree to which A. borealis may be a generalist is also important for better evaluating its invasion risk, as the establishment success of parasitoids outside their native range tends to be lower for generalists than specialists (Duncan et al. 2014).

Predicted distribution of Apocephalus borealis under future projected climate conditions
The predicted range of A. borealis changed little under the different climate change models and emission scenarios investigated (Fig. 4). Overall, future projections agreed in marginal reductions of habitat suitability at the edges of the current predicted distribution range of A. borealis and, under some scenarios (e.g., CSIRO-MK3.0 high emissions scenarios), minor increases in habitat suitability northwards. The models suggest that it is reasonable to expect the suitable habitat range of the zombie fly to remain stable until the middle of the century in spite of anticipated climate changes. However, the response of the fly to climate change will also likely depend on the availability of suitable hosts, such as Bombus bumble bees, whose range may not remain as conserved (Soroye et al. 2020).

Model evaluation and contribution of environmental variables
The AUC values for the 10 runs ranged from 0.931 to 0.935, which suggests an "excellent" (AUC > 0.90) discriminatory power of the model (Swets 1988;Araújo et al. 2005). Over the 10 cross-validation runs, the four bioclimatic variables that contributed the most to the model prediction were the maximum temperature of the warmest month (28.6.2-30.2%), isothermality (26.5-28.6%), temperature seasonality (15.4-17.6%), and annual precipitation (14.1-15.5%). The maximum temperature of the warmest month and isothermality had the highest average permutation importance (53.1-61.6% and 10.2-20.8%, respectively). The effect of each of the six studied environmental variables on the MaxEnt prediction is displayed in Fig. 4. In brief, A. borealis prefers regions with maximum temperatures of the warmest month around 20-25 °C but is rare in regions where this exceeds 35 °C. The fly likewise prefers regions with annual precipitation over 1000 mm and where precipitation of the wettest month is around 50 mm.

Is Apocephalus borealis an invasive species?
To be regarded as an invasive species, an organism must meet the following conditions: (i) become established outside of its native range, (ii) spread further from the initial site of introduction, (iii) and cause economic and/or environmental harm (Beck et al. 2008;Laštůvka and Šefrová 2020). Unfortunately, in the case of A. borealis, gaps in our knowledge of the species' biology and distribution make assessing each of these criteria difficult. The zombie fly has been reported from outside of North America on three occasions. Occurrence records from Belgium and South Korea are based on amplicons sequenced from bulk samples of managed honey bees that showed 100% similarity to the 18S ribosomal RNA sequence of A. borealis (GenBank: JF808447) (Noh et al. 2012;Ravoet et al. 2013). A BLAST search of the GenBank database conducted on 30th June 2020 by the authors however showed that there are also other phorid flies with high sequence similarity to A. borealis (≥ 95% identity), namely members of the genera Pseudacteon and Megaselia, which is not unexpected since 18S is a slow-evolving gene used in phylogenetic studies (e.g., Belinky et al. 2012). Given that 18S sequences are available for only a small share of the about 4000 known phorid fly species known to date (Brown et al. 2018), and that no specimens of the zombie fly have been collected from Asia and Europe, these reports should be treated with caution. It is more likely that the sequences originated from a common scavenger species, Megaselia scalaris in dead bees. Notably, a metatranscriptomic analysis of bee colonies in Turkey, that did not rely on single gene markers, failed to detect A. borealis in the studied samples (Tozkar et al. 2015). A paper reporting zombie flies associated with colony losses from several apiaries in Egypt was published (Khattab and El-Hosseny 2014), but based on the photographs provided, the fly in question is in fact the native Megaselia scalaris (Brown, 2019). Megaselia flies are associated with honey bees across the world, including Italy (Dutto and Ferrazzi 2014;Ricchiuti et al. 2016), Algeria (Menail et al. 2016), Cameroon (Cham et al. 2018), India (Debnath and Roy 2018), and Slovakia (Sabová 2020), and as such may account for past miss-identifications attributed to A. borealis. Nonetheless, the knowledge of the global fauna of phorids associated with Apidae remains limited (Mohammed 2018), highlighting the need for future basic research.

Reducing the invasion risk of Apocephalus borealis
Regardless of the uncertainties associated with the biology and current distribution of A. borealis, our models suggest that, if introduced outside of its native range, the zombie fly has potential to establish in vast areas of continental Europe, the Mediterranean Basin, Asia Minor, southern Africa, east Asia, Australasia, and North and South America where it could pose a threat to beekeeping and native pollinators. In the last five decades, the global apicultural sector has suffered as a result of the spread of several previously little-known invasive species, namely the Varroa mite (Varroa destructor), the small hive beetle (Aethina tumida), and the Asian hornet (Vespa velutina). All three represent major contributors to honey bee mortality (Le Conte et al. 2010;Smith et al. 2013;Oldroyd 2007). As such, much interest has recently been directed to identifying potential emerging threats to honey bee health and designing measures to prevent their introduction to other parts of the world (Oldroyd and Allsopp 2017).
The international trade in live bees is already subjected to a number of biosecurity restrictions with the aim of protecting bee health (Mutinelli 2011). The majority of international shipments of live bees from North America are destined for Australasia, South America, and to a lesser extent eastern Asia (trademap.org data for live bees, 2019), overlapping with the suitable habitat range of the zombie fly. Additional controls could be imposed on live bees originating from North America focusing on zombie flies. Although Borgmeier (1963) reported A. borealis larvae parasitising the black-widow spider Latrodectus mactans, laboratory rearing experiments were not able to verify this observation (C. Quock, personal observation), and it consequently seems that the fly's host range is restricted to hymenopterans. Since A. borealis depends on adult hymenopterans for reproduction and adult flies probably spend a considerable portion of their life locating suitable hosts, the risk of A. borealis spread could be controlled by monitoring commercially traded hymenopterans. These include honey bee packets, nucs, queen cages, and bumblebee colonies sold for pollination. The outputs of our habitat suitability model could be used to designate high-risk states and provinces, from which exports of live bees would be subjected to more rigorous checks. Bee shipments originating from high-risk A. borealis areas could be visually checked for adult flies and emerged pupae, and bee inspectors should pay close attention to bees displaying characteristically lethargic and disoriented movement and isolate them to see if A. borealis larvae emerge from them. Dead bees should be treated with ethanol or another killing agent before disposal, or disposed of by burning.
The proportions of infected bees in samples obtained from hives around the San Francisco Bay Area peaked into the summer and fall but were lower in spring. More long-duration sampling from hives in other regions of the USA and Canada is needed to determine how universal these seasonal shifts in infection prevalence might be. However, restricting the time of year when bees can be shipped to periods when parasitism is seasonally low in the region of origin might be another method of controlling A. borealis spread.
Laboratory rearing has shown that adult zombie flies seem to only search for and attack adult host during the first 3-4 days of their life. Some females reared on sugar water live to just over a week but do not attempt to seek hosts. Once infested, adult bees usually die within 7-14 days (C. Quock, personal observation). A more costly measure would thus be to quarantine shipments for 7-14 days, during which dead bees would be monitored for A. borealis emergence. Adding these preventive measures to the already existing suite of checks that transported bees are subjected to could represent a cost-effective approach to reducing the risk of future A. borealis spread.

CONCLUSIONS
The zombie fly (A. borealis) is a parasitoid of arthropods including honey bees and bumble bees that reduces the lifespan of infested insects and harbours viral and fungal pathogens (Core et al. 2012). Recently, suspected occurrences of the phorid in Belgium and South Korea uncovered during molecular screenings for bee pathogens and parasites (Noh et al. 2012;Ravoet et al. 2013) have raised concerns about the invasive potential of the species. Here, we use occurrence records of A. borealis from its native range in North America obtained from published references and the citizen science project ZomBee Watch to predict the potential distribution of the species in its native range, as well as to identify high-risk areas susceptible to further spread of the phorid and contribute to understanding its ecological requirements. Our results highlight that although at present there are no convincing records of A. borealis outside of North America, vast parts of Europe, the Mediterranean Basin, Asia Minor, southern Africa, eastern Asia, Australasia, and South America are nonetheless at high risk of invasion should the zombie fly expand beyond its native range (Fig. 5).
The global suitable habitat range is expected to remain similar at least until the middle of the twenty-first century. To mitigate the potential risk of A. borealis to world apiculture, we advise introducing simple monitoring measures of live honey bees and bumble bees exported from high risk regions of North America and highlight the need for further research of phorid parasitoids of bees.

Figure 4
Potential global distribution range of the zombie fly (Apocephalus borealis) under the future projected climate models CSIRO-MK3.0 and MIROC-H with the A1B (medium emissions) and A2 (high emissions) scenarios. Colours indicate the probability of occurrence of A. borealis (red = high, blue = low). Project, the climate modelling groups for producing and making available their model output, and two anonymous reviewers for their valuable comments. We are deeply grateful to Professor Huijie Qiao (Institute of Zoology, Chinese Academy of Sciences) for insightful comments on the manuscript. E.T. is thankful to Mattia Giacomelli (University of Bristol) for a friendly discussion regarding bioinformatics.

AUTHOR CONTRIBUTIONS
ET, SC, and CC conceived and designed research; ET conducted the modelling and analysed data; ET drafted the manuscript to which SC, CC, JH, CQ, AZ, BB, and CZ contributed; JH, AZ, CQ, and BB founded and run the ZomBee Watch project; BB is the author of the Phorid Catalog (PCAT). All authors read and approved the manuscript.

DATA AVAILABILITY
The datasets generated and analysed during the current study are available from Mendeley Data at https:// doi. org/ 10. 17632/ wvhcm wdwxk.2

CODE AVAILABILITY
Code is available from Mendeley Data at https:// doi. org/ 10. 17632/ wvhcm wdwxk.2 DECLARATIONS Ethical approval No approval of research ethics committees was required to accomplish the goals of this study because experimental work was conducted with an unregulated invertebrate species.

Conflict of interest
The authors declare no competing interests.

OPEN ACCESS
This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.