Nationwide Screening for Arthropod, Fungal, and Bacterial Honey Bee Pathogens: Utilizing Environmental DNA from Honey samples in Australia

1. Introduction

Australia is renowned for having one of the healthiest honey bee populations globally. This is largely due to its geographic isolation and strict biosecurity measures, which have kept many pests and pathogens away from the mainland [1]. However, beekeepers still need to remain vigilant for established pest and pathogenic diseases [2]. The outbreak of varroa mite in Newcastle, in June 2022 highlighted the urgent need for further research to understand better and improve our knowledge of the pathogen’s landscape affecting Australian honey bees [3].While microbial pathogens are indeed present in Australian apiaries, research on their widespread occurrence and transmission remains limited [4]. Understanding how these pathogens and parasites influence colony productivity and health requires extensive surveillance across various sites in Australia [5]. Identifying the presence of a pathogen affecting a hive can be completed by assessing hive substances, such as wax, pollen and honey. This is due to the potential transfer of pathogens and pathogen material during the production of hive substances following the acquisition of a pathogen while foraging.

In Australia, the agricultural sector heavily relies on insect pollination, particularly from the western honey bees, Apis mellifera [6] which is essential for pollinating a diverse array of crops and enhancing agricultural productivity and food security nationwide [7]. The benefits of crop pollination extend beyond agriculture exerting far-reaching impacts on the broader Australian community by facilitating pollination and positively influencing crop outcomes through their synergistic relationship with flowering plants [8]. Previous estimations indicate that the total value of paid and unpaid pollination services amounts to approximately AUD 1.2 billion annually in Australia [9]. The reliance on honeybees for optimal pollination extends to approximately two-thirds of horticultural crops in the country [10], particularly enhancing productivity and post-harvest storage qualities in fruits like apples, raspberries, and peaches [8,11,12].

These essential pollinators face significant threats from a range of bacterial, viral, and parasitic infections, which compromise their health and productivity. Environmental stressors such as exposure to harmful chemicals, inadequate nutrition and unsustainable agricultural practices also play a role [13]. Diseases such as American foulbrood (Paenibacillus larvae) and European foulbrood(Melissococcus plutonius) weaken the bees and harm brood health [14]. Infections caused by fungi, microsporidians, and arthropods contribute to colony collapse, affected by regional genetic variations and additional stressors [15,16].

Molecular methods, owing to their high sensitivity, accuracy, and capability for early pathogen detection, have been recommended by several researchers over microbial methods for detecting pathogens in environmental samples [17,18]. Polymerase chain reaction-based assays have been developed for the detection of honey bee pathogens from bees and other hive materials [19,20]. Conventional PCR is the most widely used method for detecting bacteria (Paenibacillus larvae, Melissococcus plutonius), fungi (Nosema spp., Ascosphaera apis), and arthropods (Aethina tumida and Galleria mellonella) from honey and other hive samples [21,22,23,24]. PCR is highly regarded for its consistent, specific, and precise identification of a multitude of pathogens, making it an essential tool for routine pathogen screening in surveillance operations [25]. Traditionally, adult bees have been the primary source of pathogen detection [26,27], the evolution of environmental DNA/RNA-based detection techniques presents a superior alternative to direct sampling from hosts [28]. This advancement has prompted numerous research initiatives employing PCR-based assays to detect bee pathogens and parasites directly from honey samples, offering a less invasive yet effective surveillance approach [29,30].

Monitoring disease in adult bees is challenging due to the resource-intensive requirements and costs associated with collecting data from individual bees. This method provides only a snapshot of disease prevalence at specific times, limiting the ability to track changes over time [31]. In contrast, hive materials such as honey provide a valuable source of environmental DNA (eDNA) for the detection of invasive organisms.

Honey eDNA has thus emerged as a practical method for monitoring honey bee pathogens and parasites [32,33]. Identifying the presence of a pathogen affecting a hive can be completed by assessing hive substances, such as wax, pollen, and honey. This is due to the potential transfer of pathogens and pathogen material during the production of hive substances following the acquisition of a pathogen while foraging. The genetic material left by organisms in the environment, known as environmental DNA (eDNA), is a persistent biomolecular marker that can be collected, extracted, and analyzed from various substrates, making it a powerful tool for detecting and monitoring both microbial and macrobial communities effectively [34,35]. Honey’s stable nature, characterized by low water content and acidic pH, preserves remnants of microorganisms, offering insight into historical records of colony pathogens[36]. Recognizing honey as a reservoir of exogenous DNA underscores the importance of monitoring and management of bee health, with the detection of pathogens from honey serving to prevent the spread of diseases among colonies.

In our study, we collected honey samples from diverse botanical and geographical regions spanning different states in Australia. These samples were utilized to isolate eDNA using the bead-beating-silica DNA extraction method and identified pathogens and parasites known to impact honey bee populations in Australia. We also examined the co-occurrence pattern of these pathogens in honey samples.

2. Materials and Methods

2.1. Acquisition of Commercial Honey Samples

Honey samples were gathered between 2022 and 2024 from diverse regions across Australia, including trade markets and directly from beekeepers, to ensure a comprehensive representation of honey sources. The sample comprised both poly-floral and mono-floral varieties, were obtained from different states: Victoria (27), New South Wales (29), Queensland (24), Northern Territory (1), Western Australia (22), South Australia (20* including 6 samples from kangaroo island), Tasmania (12)[37].

2.2. eDNA Extraction from Honey Samples

DNA was extracted from honey following an in-house extraction protocol, modified from Waiblinger H-U, et al., [38] with adjustments, incorporating both pre-treatment and a post-treatment phase as detailed by Soares et al., [39] and Rathinasamy et al., [40].

For the pre-treatment phase, 50 g of honey was evenly partitioned into four sterile 50 mL conical tubes, each containing 12.5 g of honey. Subsequently, 40 mL of ultrapure water was added to each tube and vortexed until completely homogenized. The samples underwent incubation at 40°C for 10 mins in a water bath with agitation, then centrifugation at 4700× g for 35 minutes. The resulting supernatant was discarded, and the pellets were resuspended in 5 mL of ultrapure water and combined with the same honey sample into a single conical tube. This suspension was then subjected to an additional centrifugation at 4700× g for 30 minutes. After centrifugation, the supernatant was removed, and the pellet was resuspended in approximately 500 µL of ultrapure water. Subsequently, the suspension was transferred to a 2 mL tube containing seven glass beads (approximately 5 mm in size) and vortexed for 2 minutes. The glass beads were removed, and the mixture underwent centrifugation at 11,000× g at 4°C for 10 minutes. The resulting pellet served as the starting material for DNA extraction.

To each pre-treated sample pellet, 860 µL of TNE (10 mM Tris-HCL, 150 mM NaCl, 2mM ethylenediaminetetraacetic acid, 1% (w/v) sodium dodecyl sulfate, pH 7.5) extraction buffer was added, which was pre-heated at 60°C. Additionally, 100 µL of 5 M guanidine hydrochloride (w/v) and 40 µL of proteinase K solution (20 mg/ml^-1) were added and subsequently vortexed. Following a 3-hour incubation at 60°C in a Thermomixer Comfort (Eppendorf AG, Hamburg, Germany) with agitation at 900 rpm, the suspension underwent centrifugation for 15 minutes at 17,000×g at 4°C. The supernatant was collected for DNA purification.

2.3. eDNA Purification

The resulting extracted DNA was mixed with 2 volumes of 6M sodium iodide (NaI) and 100 µL of 100 mg/mL silica dioxide (SiO₂), and binding occurred with gentle agitation on a rocker for 30 minutes. Subsequently, centrifugation was performed (10 minutes at 5000×g at 4° C), with the supernatant discarded and 500 µL of silica wash buffer (50% (v/v) ethanol, 10 mM Tris-HCl, 100 mM NaCl, 1 mM EDTA, pH 8) added to the pellet and resuspended via vortexing. The solution was centrifugated for 1 minute at 4700×g at 4°C, with this wash step repeated a further two times. Care was taken to remove all the supernatants. To elute the eDNA from the silica matrix, 50 µL of elution buffer (10 mM Tris-HCl, pH 8) was added and incubated at 70°C for 5 minutes. Following incubation, the sample was centrifugated for 5 minutes at 16,000×g, with the supernatant collected and stored at -20°C until further analysis. The concentration of DNA was measured using a Nanodrop Eight spectrophotometer and visualized by electrophoresis in a 2% (w/v) agarose gel in 1X TBE buffer (0.13 M Tris, 45 Mm Boric acid, 2.5 mM EDTA pH 7.6) containing 0.2 µg/mL ethidium bromide (Merck Life Science Pty Ltd., Melbourne, Australia).

2.4. PCR Analysis

Gene targets and primer sequences were selected based on validated and published studies, as outlined in Table 1. We assessed the presence of pathogens using both singleplex and multiplex PCR, with PCR conditions optimized for each pathogen according to their annealing temperature and amplicon size (Table 2). Positive DNA controls for each pathogen were generated by chemical synthesizing of each gene (Table 1) from Integrated DNA Technologies (IDT) and subsequently cloning them into the TOPO vector according to manufacturer instructions. Sensitivity was determined by serial dilution of synthetic plasmid from 10^-1 ng/µl to 10^-9 ng/µL concentration, and specificity was evaluated using genomic DNA of P. larvae, M. plutonius, A. apis, N. apis, N. ceranae, A. tumida, and G. mellonella.

PCR reactions were performed in a total volume of 25 µl, containing 1x GoTaq^® Green Master Mix (Promega), 25ng of template DNA, and species-specific primer concentration. The primer concentrations were optimized as follows: 0.5 µM for A. mellifera, N. apis, N. ceranae, and G. mellonella; 0.4 µM for P. larvae and M. plutonius; 0.3 µM for A. tumida; and 0.1 µM for A. apis. The thermal cycling conditions were performed as described in Table 2. A negative control, devoid of template DNA, was included in each assay. Subsequently, the amplified DNA fragments were subjected to electrophoresis in a 2% (w/v) agarose gel in 1X TBE buffer containing 0.2µg/mL ethidium bromide. Three confirmed positive samples for each pathogen were sequenced to verify the correct target amplification.

2.5. Statistical Analysis

The diversity of pathogens detected in each honey sample was quantified and visualized using the ggplot2 package in R [41]. To assess the patterns of pathogen co-occurrence across the samples, we utilized the R package cooccur [42], which employs a probabilistic model to evaluate species co-occurrence [43]. This model includes combinatorial methods that calculate the probability of the observed frequency of co-occurrence being significantly greater (p_gt < 0.05) than expected, indicating a positive association (+ve) and with p_gt > 0.05 indicating negative association (-ve), or independent co-occurrence (p_lt > 0.05), which indicates no evidence of pathogen interaction and suggests no significant association. The effect size, which ranges from -0.02 to 0.08, indicates the strength of co-occurrence (+ve/-ve) between two pathogens present in a sample and is assessed using the probability matrix (p_gt and p_lt). In total, 21 species pairs were analyzed to determine the frequency of their co-occurrence across the samples and the nature of the association (+ve or -ve) between each pathogen pair.

Table 1. PCR primers used in this study to amplify eDNA extracted from honey samples.

Target species	Primer name1	Accession No.	Primer sequence (5′-3′)	Amplified region	Product size (bp)	Reference
Singleplex PCR
Apis mellifera	AM Forward AM Reverse		GGCAGAATAAGTGCATTG TTAATATGAATTAAGTGGGG	suppl	C 85, M 1392	[44]
Nosema apis	Nose_apis_chen_F Nose_apis_chen_R	U97150.1	CCATTGCCGGATAAGAGAGT CCACCAAAAACTCCCAAGAG	SSUrRNA	269	[45]
Nosema ceranae	Nose_cera_chen_F Nose_cera_chen_R	DQ486027.1	CGGATAAAAGAGTCCGTTACC TGAGCAGGGTTCTAGGGAT	SSUrRNA	250	[45]
Aethina tumida	Atum-3F Atum-3R	MF943248.1	CCCATTTCCATTATGTWYTATCTATAGG CTATTTAAAGTYAATCCTGTAATTAATGG	COI	97	[46]
Galleria mellonella	GallMelCox1-F GallMelCox1-R	KT750964.1	TGAACTTGGTAATCCTGGTTCT TATTATTAAGTCGGGGGAAAGC	COI	182	[46]
Multiplex PCR
Paenibacillus larvae Melissococcus plutonius Ascosphaera apis	Han233PaeLarv16S_F Han233PaeLarv16S_R Mp_Arai187_F Mp_Arai187_R AscosFORa AscosREVa	NZCP019687.1 AB778538.1 U68313.1	GTGTTTCCTTCGGGAGACG CTCTAGGTCGGCTACGCATC TGGTAGCTTAGGCGGAAAAC TGGAGCGATTAGAGTCGTTAGA TGTGTCTGTGCGGCTAGGTG GCTAGCCAGGGGGGAACTAA	16S rRNA NapA 18S rRNA	233 187 136	[47]

¹ The internal name of the forward and reverse primers. ² The size of the amplified fragment may vary depending on the mitochondrial lineage C which is highly frequent in A.m. lingustica and M in A.m. mellifera (C or M).

Table 1. PCR primers used in this study to amplify eDNA extracted from honey samples.

Target species	Primer name1	Accession No.	Primer sequence (5′-3′)	Amplified region	Product size (bp)	Reference
Singleplex PCR
Apis mellifera	AM Forward AM Reverse		GGCAGAATAAGTGCATTG TTAATATGAATTAAGTGGGG	suppl	C 85, M 1392	[44]
Nosema apis	Nose_apis_chen_F Nose_apis_chen_R	U97150.1	CCATTGCCGGATAAGAGAGT CCACCAAAAACTCCCAAGAG	SSUrRNA	269	[45]
Nosema ceranae	Nose_cera_chen_F Nose_cera_chen_R	DQ486027.1	CGGATAAAAGAGTCCGTTACC TGAGCAGGGTTCTAGGGAT	SSUrRNA	250	[45]
Aethina tumida	Atum-3F Atum-3R	MF943248.1	CCCATTTCCATTATGTWYTATCTATAGG CTATTTAAAGTYAATCCTGTAATTAATGG	COI	97	[46]
Galleria mellonella	GallMelCox1-F GallMelCox1-R	KT750964.1	TGAACTTGGTAATCCTGGTTCT TATTATTAAGTCGGGGGAAAGC	COI	182	[46]
Multiplex PCR
Paenibacillus larvae Melissococcus plutonius Ascosphaera apis	Han233PaeLarv16S_F Han233PaeLarv16S_R Mp_Arai187_F Mp_Arai187_R AscosFORa AscosREVa	NZCP019687.1 AB778538.1 U68313.1	GTGTTTCCTTCGGGAGACG CTCTAGGTCGGCTACGCATC TGGTAGCTTAGGCGGAAAAC TGGAGCGATTAGAGTCGTTAGA TGTGTCTGTGCGGCTAGGTG GCTAGCCAGGGGGGAACTAA	16S rRNA NapA 18S rRNA	233 187 136	[47]

¹ The internal name of the forward and reverse primers. ² The size of the amplified fragment may vary depending on the mitochondrial lineage C which is highly frequent in A.m. lingustica and M in A.m. mellifera (C or M).

Table 2. Optimized PCR conditions for detection of pathogens.

Target Species	Steps	Optimized Conditions	Time	Cycle
Multiplex PCR	Initial Denaturation Annealing Extension Final Extension	95°C 95°C 63°C 72°C 72°C	2 min 1 min 1 min 1 min 5 min	35x
P. larvae M. plutonius A. apis
Single plex PCR
N. apis N. ceranae	Initial Denaturation Annealing Extension Final Extension	95°C 94°C 58.6°C 68°C 72°C	2 min 15 sec 30 sec 1 min 7 min	35x
A. tumida	Initial Denaturation Annealing Extension Final Extension	95°C 98°C 54°C 72°C 72°C	3 min 20 sec 30 sec 1 min 7 min	35x
G. mellonella	Initial Denaturation Annealing Extension Final Extension	95°C 98°C 61°C 72°C 72°C	3 min 1 min 1 min 1 min 1 min	35x

Table 2. Optimized PCR conditions for detection of pathogens.

Target Species	Steps	Optimized Conditions	Time	Cycle
Multiplex PCR	Initial Denaturation Annealing Extension Final Extension	95°C 95°C 63°C 72°C 72°C	2 min 1 min 1 min 1 min 5 min	35x
P. larvae M. plutonius A. apis
Single plex PCR
N. apis N. ceranae	Initial Denaturation Annealing Extension Final Extension	95°C 94°C 58.6°C 68°C 72°C	2 min 15 sec 30 sec 1 min 7 min	35x
A. tumida	Initial Denaturation Annealing Extension Final Extension	95°C 98°C 54°C 72°C 72°C	3 min 20 sec 30 sec 1 min 7 min	35x
G. mellonella	Initial Denaturation Annealing Extension Final Extension	95°C 98°C 61°C 72°C 72°C	3 min 1 min 1 min 1 min 1 min	35x

3. Results

3.1. Assessment of Extracted DNA

Before attempting to analyze for the presence or absence of honey bee pathogens, we optimized PCR for the amplification of various pathogens using either single- or multiplex PCR. The revalidation of each primer set was completed to ensure high sensitivity and specificity. Sensitivity is determined at the lowest detectable limit (10^-5 ng/µL), with no cross-reactivity observed when testing genomic DNA from P. larvae, M. plutonius, A. apis, N. apis, N. ceranae, A. tumida, and G. mellonella (Supplementary Fig. S1).

To verify the successful extraction and purification of eDNA from all honey samples, the presence of the Apis mellifera mtDNA was determined by PCR. The amplification resulted in fragments of 85 base pair (bp), specific to A. mellifera C lineage prevalent in A. m. lingustica and 139 bp specific to the M lineage, characteristic of A. m. mellifera (Figure 1). Successful PCR amplification was achieved for all tested honey samples (Supplementary Figure S2), indicating the successful extraction and purification of eDNA from the honey.

Every honey sample underwent PCR analysis to detect various common pathogens and pests that affect honey bees (Supplementary Figure S3). Specifically, samples 107 and 108 were positive for N. apis, while all samples except 104 tested positive for N. ceranae (Figure 1A). In the multiplex PCR targeting bacterial pathogens P. larvae, M. plutonius, and fungal pathogen A. apis, sample 107 amplified all three pathogens, whereas sample 108 tested positive for P. larvae and M. plutonius (Figure 1A). No amplification was detected in samples 104, 105, 106, 109, and 110 (Figure 1A). Furthermore, samples 125-131 indicated the presence of the arthropod pests A. tumida and G. mellonella (Figure 1B). All samples tested positive for A. tumida except samples 126 and 130, which tested positive for G. mellonella (Figure 1B).

3.2. Prevalence Pattern Across Different Australian States

Analysis of 135 honey samples collected from various locations in Australia revealed N. ceranae emerged as the most prevalent pathogen, present in 57% of the samples. This was followed by A. tumida (40%), G. mellonella (37%), P. larvae (21%), N. apis (19%), and M. plutonius (18%). A. apis was detected in a smaller proportion of the samples, with a prevalence of 5%. Additionally, 25 samples (19%) tested negative for all pathogens analysed (Figure 2). The most recent national survey of honey bee pathogens and pests, last conducted in 2014, reveals a significant shift in pathogen levels and their distribution throughout Australia over the following seven years. The detailed geographical distribution of bacterial, fungal pathogens, and arthropod parasites throughout Australia over the two years is summarized in Table 3.

Regional variations in pathogen prevalence were evident across different Australian states. Specifically, American foulbrood (AFB) and European foulbrood (EFB) were identified in Victoria (VIC), New South Wales (NSW), Queensland (QLD), and Tasmania (TAS) (Table 3). Notably, Victoria exhibited a high prevalence, with 33% of the sample testing positive for AFB and 37% for EFB. In NSW and QLD, there was no significant difference in the prevalence of AFB (28% and 21%, respectively) and EFB (14% and 13%, respectively). Tasmania had a prevalence of 25% for both AFB and EFB. In contrast, Western Australia only harbours AFB (14%) and is free from EFB, while South Australia exclusively hosts EFB (29%), with AFB absent. Similarly, EFB was absent in Western Australia (WA) and Kangaroo Island (KI) (Table 3).

The primary fungal disease impacting honey bees, A. apis (commonly referred to as chalkbrood), was detected at low levels in Victoria (4%) and New South Wales (7%), with no significant presence in Queensland, Northern Territory, Western Australia, and South Australia (Table 4). In contrast, Tasmania exhibited a relatively high incidence of A. apis at 33% (Table 4).

Nosema ceranae was frequently found across all states, with a particularly high prevalence in VIC (70%), followed by SA (64%), TAS (50%) and WA (45%). Remarkably, even a single sample from NT tested positive for N. ceranae (Table 4). There was no significant variance in prevalence between NSW (59%) and QLD (58%). N. apis was identified in TAS (58%), WA (23%), VIC (19%), SA (14%), and NSW (10%), but it was not detected in QLD or NT (Table 4). Additionally, 16% of the samples had mixed infections with both N. ceranae and N. apis (Table 4).

The invertebrate pests of honey bees, A. tumida and G. mellonella, were present in honey samples from all states except KI (Table 5). A. tumida showed a high prevalence in QLD (71%). SA (57%), VIC (56%), TAS (33%), and NSW (31%), while WA had a significantly lower prevalence at 5%. On the other hand, G. mellonella was most prevalent in TAS (84%) followed by SA (57%), WA (45%), VIC (44%), QLD (25%) and NSW (10%). Notably, in SA, 57% of the samples tested positive for both A. tumida and G. mellonella (Table 5).

3.3. Pathogen Prevalence on Kangaroo Island

Kangaroo Island is a distinct region due to its island status and strict movement restrictions, which help preserve the pure bred honey bees, preventing crossbreeding. Six samples were collected from the commercial beekeepers on KI and subjected to testing for bacterial, fungal pathogens, and arthropod parasites. The results were significant, as all samples tested negative for the three brood diseases (P. larvae, M. plutonius, and A. apis), as well as the pests A. tumida and G. mellonella. The only honey bee pathogens recorded on KI are the fungal pathogens N. apis and N. ceranae, with 50% and 17% of samples testing positive, respectively (Table 4).

3.4. Trends in Concurrent Infections

This study uncovered a diverse array of honey bee pathogens and pests present in individual honey samples, highlighting a significant incidence of co-infection within the bee population. The analysis revealed that 30% (40/135) of the samples contained a single type of pathogen, while 20% (27/135) displayed two distinct pathogens. Co-infections involving 3 and 4 pathogens were identified in 13% (18/135) and 10% (13/135) of the samples, respectively. Furthermore, 6% (8/135) contained 5 different types of pathogens, with less than 4% (6/135) having 6 types of infection. Notably, none of the samples tested positive for all seven types of pathogens analyzed (Figure 4A).

The probabilistic co-occurrence model revealed several significant positive associations (Figure 4B). Among these, positive co-occurrences were observed between P. larvae and M. plutonius (p_gt = 0.00001, p_lt =1, effect size = 0.06666 ). Notably, M. plutonius exhibited strong to moderate positive co-occurrence with all studied pathogens, including moderate association between A. apis (p_gt = 0.00198, p_lt = 0.99989, effect size = 0.02814), N. apis (p_gt = 0.02280, p_lt = 0.99389, effect size = 0.03111), N. ceranae (p_gt = 0.02776, p_lt = 0.99171, effect size = 0.03481), and a strong positive association with A. tumida (p_gt = 0.00082, p_lt = 0.99985, effect size = 0.05481), and G. mellonella (p_gt = 0.00017, p_lt = 0.99997, effect size = 0.06148). Additionally, A. apis displayed weak positive co-occurrence with N. apis (p_gt = 0.02949, p_lt = 0.99639, effect size = 0.01925) and G. mellonella (p_gt = 0.00931, p_lt = 0.99938, effect size = 0.02592), while N. apis showed moderate positive association with N. ceranae (p_gt = 0.00752, p_lt = 0.99812, effect size = 0.04444) and G. mellonella (p_gt = 0.01889, p_lt = 0.99403, effect size = 0.03851). N. ceranae had a strong positive association with A. tumida (p_gt = 0.00120, p_lt = 0.99966, effect size = 0.06666) and a moderate positive association with G. mellonella (p_gt = 0.00407, p_lt = 0.99873, effect size = 0.05777). A strong positive co-occurrence was also observed between the two arthropod pests (p_gt = 0.00058, p_lt = 0.99985, effect size = 0.06962).

Among the 21 pathogen pairs examined, 13 pairs exhibited positive associations, ranging from weak to strong, and tended to co-occur more frequently than expected (Supplementary Table S2). In contrast, eight pathogen pairs showed no statistically significant positive or negative association, suggesting no evidence of pathogen interaction. For instance, the interaction between P. larvae and N. ceranae is characterized by a weak negative association (p_gt = 0.67507, p_lt = 0.48867, effect size = -0.00444), indicating that the two pathogens are likely independent of each other. Similarly, N. apis and A. tumida exhibited a weak negative association (p_gt = 0.28640, p_lt = 0.84380, effect size = -0.01333). Additionally, P. larvae showed no positive or negative association with A. apis (p_gt = 0.15581, p_lt = 0.96632, effect size = 0.01111), N. apis (p_gt = 0.06581, p_lt = 0.97716, effect size = 0.02518), A. tumida (p_gt = 0.28489, p_lt = 0.84055, effect size = 0.01333) or G. mellonella ( p_gt = 0.15120, p_lt = 0.92831, effect size = 0.02074). Similarly, A. apis did not exhibit significant co-occurrence with N. ceranae ( p_gt = 0.67771, p_lt =0.62275, effect size = 0.00074) or with A. tumida ( p_gt = 0.28532, p_lt = 0.90959, effect size = 0.00888) showing that these pathogen pairs occur independently of each other with no strong evidence of interaction.

4. Discussion

Honey serves as an ideal source of environmental DNA because foraging honey bees collect and transfer environmental microorganisms and contaminants to the hive, depositing them in the honeycomb and honey [30,33]. This eDNA reveals the presence of various organisms, including bee pathogens and parasites, making honey a valuable non-invasive tool for monitoring colony health and detecting invasive species [21,48]. Many studies have developed PCR-based diagnostic assays to detect key honey bee pathogens and pests in honey samples, including Paenibacillus larvae, Melissococcus plutonius [49,50], and Ascosphaera apis [47]. Giersch et al., identified Nosema apis and Nosema ceranae in adult bees and honey samples from Australia (2007–2008) using PCR, Restriction fragment length polymorphism (RFLP), and microscopy [51]. Later, Granato et al., employed PCR to detect N. apis and N. ceranae spores in honey, providing a reliable method to track infections in bee colonies, even without visible clinical symptoms [52]. Expanding on this approach, Ribani et al. developed a DNA-based assay to detect traces of Aethina tumida and Galleria mellonella in Italian honey samples [46] supporting honey as a source of eDNA to analyse and detect hive-associated organisms. Further investigation by Ribani et al., and Aditya et al., [18,21,46] used eDNA in honey to detect specific pathogens and found it to be an effective method for the global monitoring of pathogens.

Our data highlights the presence of key honey bee pathogens in honey samples collected across Australia, with Victoria, New South Wales, and Tasmania showing the highest prevalence of pathogens and parasites compared to other states. The most widely distributed and frequently detected pathogens include microsporidian N. ceranae, the arthropod pest A. tumida, and G. mellonella. The bacterium P. larvae, M. plutonius, and microsporidian N. apis were detected at low prevalence, with A. apis detected at low frequency. These findings provide novel trends in pathogen prevalence and distribution and identify key disease hotspots in Australia. The presence of these pathogens in honey samples may indicate residual DNA from past infections or ongoing infections [53].

Nosemosis, caused by the microsporidian fungal pathogens N. apis and N. ceranae, is a globally distributed honey bee disease [54]. N. ceranae was initially believed to be restricted to Asian honey bees (Apis ceranae) when first detected in Beijing, China, in the 1990s [55]. However, A. mellifera was previously believed to be susceptible only to N. apis. This view changed when N. ceranae was identified in A. mellifera during a study in Spain, using isolates from various regions worldwide [56,57]. The host shift of N. ceranae from A. ceranae to A. mellifera has occurred on a global scale, raising questions about the factors driving this transition, but likely due to the movement of infected honey bees through commercial trade [58,59,60,61]. The high prevalence of N. ceranae compared to N. apis in Australian honey samples is part of a broader trend observed in other countries. This result aligns with surveys that report a high prevalence of N. ceranae in Canada [62], in the United States [63], Florida [64], Mexico [65], Northern and Central Italy [66,67], Scotland [68], New Zealand [69,70,71] and Indonesia [18]. N. ceranae was isolated from diseased A. mellifera in Taiwan, where it is considered a non-native species [72]. However, studies in Germany [73] and Sweden [74] reported a higher occurrence of N. apis compared to N. ceranae, likely due to climatic factors. Experimental Studies indicate N. ceranae exhibits reduced germination and infectivity at low temperatures, while N. apis remain viable even under freezing conditions [75]. This differential temperature tolerance likely contributes to the higher prevalence of N. apis in colder regions, whereas N. ceranae flourishes in warmer areas [76].

The analysis conducted seven years ago on 151 adult bee samples in Australia revealed that 139 of the samples tested positive for either N. apis or N. ceranae, with an equal proportion of both pathogens. Interestingly, there is a gradual decline in the prevalence of N. apis among the samples tested in this study, indicating a shift in pathogen dynamics over time. Giersch et al., [51] detected N. apis in bee samples across various Australian states, while N. ceranae was found in all states except Western Australia and Tasmania. In contrast, our study observed N. ceranae in all sampled states, but N. apis was not detected in Queensland. This higher prevalence of N. ceranae compared to N. apis in this study may be attributed to Australia’s warmer temperatures, which favor N. ceranae, making it a more competitive pathogen over time. The occurrence of co-infection was more frequent in tested honey samples for N. apis and N. ceranae than single-species infections, with N. ceranae exhibiting dominance over N. apis. Similar findings of mixed infections have been reported in Europe, North America, and Asia [77]. The ability of N. ceranae to thrive in warmer conditions, its rapid reproduction rate, and significantly higher mortality in honey bees suggest its competitive advantage over N. apis [60]. N. ceranae is widespread across all states in Australia, which have diverse climate conditions ranging from tropical in the north to temperate in the south. Further investigation is required to substantiate the correlation between climate and the prevalence of Nosema spp.

Chalkbrood was detected in 5% of the analyzed honey samples, restricted to Victoria and New South Wales. Although A. apis was first reported in Queensland, no traces were found in honey samples from this region. Robert et al., [2] identified chalkbrood mummies in 66% of hives across Australia, primarily along the eastern coast, with a single case recorded in Western Australia. However, our study did not identify A. apis in honey samples from this state. Globally, chalkbrood is widely distributed and reported in key beekeeping regions, including Central America, North America, Mexico, Chile, Japan, China, Turkey, Africa, and the Philippines [78,79]. The spores of A. apis, N. apis, and N. ceranae remain viable for over 15 years in infected apiaries, stored honey, hive materials, and beekeeping equipment, serving as a persistent source of infection [80,81]. Determining whether the detected genetic material or isolated DNA comes from viable spores or non-viable residue or is transmitted by exchanging contaminated honey and hive components between colonies is important.

American Foulbrood and European Foulbrood are serious bacterial diseases caused by Paenibacillus larvae and Melissococcus plutonius, are widespread, highly contagious and can lead to significant losses in honey bee colonies [82]. Both AFB and EFB are globally distributed honey bee diseases. AFB is present in all beekeeping countries and is classified by the World Organization for Animal Health (WOAH) as a highly dangerous infectious disease in animals. First identified in North America [83,84] and its causative agent, P. larvae, produces resilient, long-lived spores in large quantities, which are essential for infection. Once the colony is infected, recovery is not possible [85]. Whereas EFB has long been recognized in Europe and North America, more recently, it has also been detected in Africa, South America, India, Japan, and Australia [86].

The prevalence of AFB and EFB varies across regions due to differences in sampling methods, colony health, and environmental factors. In Australia, Robert et al., [2] reported AFB at 9% and EFB at 5% in 2015. However, this study found a higher nationwide prevalence, with 21% for AFB and 18% for EFB. This difference is likely due to the presence of spores in honey samples, which may persist for an extended period, compared to samples collected directly from hives exhibiting clinical symptoms of AFB and EFB. These results remain consistent over time, with WA, NT, and KI continuing to be free from EFB. In Indonesia, P. larvae was not detected, whereas M. plutonius was identified for the first time from honey samples collected from Java [18]. Ribani et al., [21] reported a high prevalence of EFB (87%) and AFB (49%) in Italian honey samples processed between 2004 and 2018. Similarly, P. larvae, spores were detected in 40% of the analyzed samples in Italy, indicating a long-term persistence of the disease [87]. In the U.S, EFB (19.2%) was more common in symptomatic colonies than AFB (8.6%), despite the resilience and ease of spread of P. larvae spores. The higher prevalence of EFB is likely due to its tendency to infect already stressed colonies [88]. Given the long-term viability of P. larvae spores (up to 15 years), the risk of disease spread through honey exchange remains a major concern, necessitating ongoing surveillance and strict biosecurity measures [85].

Arthropod pests A. tumida and G. mellonella are among the most damaging threats to honey bees, particularly when they spread beyond their native habitats [89]. A. tumida, originally native to Sub-Saharan Africa, and G. mellonella, once limited to A. ceranae, have also expanded their presence worldwide [90]. The larvae stage of both pests harm the colony, with A. tumida causing hive destruction and honey fermentation [91], whereas wax moth larvae consume wax combs, honey, and the cast skin of bee larvae [92]. A. tumida and G. mellonella are widely distributed in all continents except Antarctica [93,94,95]. The small Hive Beetle (SHB) was first detected in the U.S. [96] and later in Australia [97]; thereafter, it spread to all parts of the world. In South Africa, 69% of apiaries tested positive for SHB and wax moth infestation between 2010 and 2011 [98]. In Nigeria, apiaries are experiencing a concerning infestation rate of 21% for SHB and 5% for wax moths [99]. The comparison of SHB and wax moth prevalence in adult bee samples from 2015-2017 [2] with current honey samples in Australia reveals a significant trend. In the earlier study, SHB was detected to have a high incidence in Queensland but was absent in South Australia, Western Australia, Tasmania, and Kangaroo Island. At the same time, wax moth was present across all states. The honey eDNA analysis in our study detected both SHB and wax moths in all states except Kangaroo Island, with SHB still being more prevalent in Queensland. This suggests that the range of SHB may be expanding, which could have significant implications for pest management, while the distribution of wax moths remains extensive [100,101].

The investigation of pathogen co-infection in bee hives is limited, with few studies addressing the co-occurrence pattern of various honey bee pathogens and pests. Our data reveals a strong to modrate positive co-occurrence of M. plutonius with all the tested pests and pathogens, consistent with findings by Deutsch, Kaitlin R., et al., [64] which linked M. plutonius to A. apis and N. ceranae. Aglagane et al.,[102] reported higher prevalence levels of N. ceranae and M. plutonius in migratory beehives compared to stationary colonies. M. plutonius remains viable in honey, beeswax, and pollen for years without any clinical symptoms [103,104] and often kills honey bee larvae, leading to a decreased adult honey bee population and weakening the colony strength and productivity, exposing the bees to other pest and pathogen attacks [23,105]. A strong positive association between both bacterial diseases has been reported by Sturtevant et al., [106], indicating the presence of both bacterial species in the honey bee colonies. One possible explanation for this strong positive co-occurrence is that the presence of AFB can affect overall health and immune defense, creating an environment where EFB can easily take hold [106,107].

A study in 32 Kenyan apiaries revealed a significant positive correlation between the abundance of A. tumida and N. ceranae using the Spearman correlation coefficient (P = 0.01507) [108] . Our study, using the probabilistic model (R package cooccur), also demonstrated a positive correlation between A. tumida and N. ceranae with a p_gt = 0.00120, which is below the established significant threshold of 0.05, and an effect size = 0.06666, supporting the strength of association. This suggests a synergistic relationship; as the number of A. tumida increases, the prevalence of Nosema infection also tends to rise. This evidence suggests the potential role of A. tumida in the biological transmission of N. ceranae, as hive beetles can carry spores from infected bee feaces [109]. Extensive research has been conducted on the co-infection of Nosema apis and Nosema ceranae in honey bee colonies, demonstrating correlations with seasonal variation and asymmetrical coexistence [45,110,111]. Our study corroborates these findings, confirming a strong positive association between the two Nosema species.

The trends in pathogen co-occurrence observed in this study highlight notable patterns among different pests and pathogens. However, these patterns do not confirm that these pathogens directly cause infections or facilitate each other’s spread, indicating a knowledge gap in understanding the co-occurrence of honey bee pathogens and pests. This data provides an interesting foundation for further exploration in future studies.

5. Conclusions

For the first time, a detailed pathogen prevalence survey using honey from different states of Australia revealed that N. ceranae is the most prevalent pathogen, followed by A. tumida, G. mellonella, P. larvae, N. apis, M. plutonius, and A. apis. In particular, detecting nosemosis in Kangaroo Island indicates the need to address emerging threats to honey bee populations. Our data on the current prevalence of arthropod, fungal, and bacterial pathogens in honey bees can contribute to safeguarding Australia’s beekeeping industry.