Using whole-genome sequencing to assess the diversity of Paenibacillus larvae within an outbreak and a beekeeping operation.

Data Summary

The sequencing data have been deposited in the NCBI Short Sequence Archive (SRA) database under the BioProject accession number PRJNA725965. A complete list of the genomes analysed in this study is available in Table S1 at Figshare: https://doi.org/10.6084/m9.figshare.16573895.v1 [1].

In this work, we describe an outbreak of American foulbrood (AFB) in Slovenia in 2019, which was caused by two ERIC II-ST11 outbreak clones. Genomic data were interpreted in combination with epidemiological data. The median distance between the isolates from a single beekeeping operation was 10 whole-genome multilocus sequence typing allele differences (AD; range=1–22), similar to the diversity observed in individual apiaries and colonies. The maximum pairwise AD for the two outbreak clusters was 33 and 44; the outbreak-related isolates differed from the non-outbreak isolates in a minimum of 37 AD. The present study provides insight into genetic diversity of at different levels, providing information for future whole genome sequencing-based studies of AFB and emphasizing the importance of the associated epidemiological data when delineating outbreaks based on genetic thresholds.

Introduction

American foulbrood (AFB) is a devastating disease of honeybees (Apis mellifera) and other species of the genus Apis caused by [2]. According to Council Directive 92/65/EEC [3], AFB is a notifiable disease in the European Union and requires strict control and eradication measures, resulting in considerable economic losses to the beekeeping sector. In the past two decades, conventional PCR based on the enterobacterial repetitive intergenic consensus (ERIC) sequences [4] has become the main genotyping method for , classifying the species into five ERIC types (ERIC I–V) [5-7]. Only ERIC I and II are epidemiologically relevant types and are distributed globally [2, 8]. ERIC I exhibits a slow-killing phenotype (after comb cell capping) and is thus more virulent at the colony level, whereas ERIC II shows a fast-killing phenotype with increased virulence at the larval level [6, 7, 9].

Whole-genome sequencing (WGS) has revolutionized the genomic epidemiology of pathogenic microorganisms by providing unprecedented discriminatory power and detailed insights into gene content and diversity. Three widely used WGS-based approaches are core-genome multilocus sequence typing (cgMLST), whole-genome MLST (wgMLST) and SNP-based typing [10-12]. Whereas cg/wgMLST are both allele-based approaches and require the construction of a cg/wgMLST scheme, SNP typing is based on read mapping against a closely related reference genome. Recently, stable cg/wgMLST schemes for have been developed, aiming to harmonize WGS-based typing of [13, 14]. We have previously shown that all three approaches provide generally consistent results in terms of delineating outbreak clusters [14].

WGS-based outbreak investigation relies on the identification of clusters of closely related (clonal) isolates. Only a limited number of studies have applied WGS for outbreak investigation of AFB and aimed to establish a threshold to distinguish between outbreak-related and non-related isolates [13-16]. Nonetheless, establishing an optimal threshold with maximum sensitivity to include the outbreak-related isolates and maximum specificity to exclude the non-outbreak isolates is a non-trivial task that requires integration of genomic and epidemiological data [11, 12, 17]. Thus, the suggested thresholds serve as orientation, but each investigated case should be interpreted on a case-by-case basis and in conjunction with the available epidemiological data. Further research is needed to re-evaluate the proposed thresholds for .

The aim of the present study was to evaluate the intra- and inter-outbreak diversity of isolates with particular emphasis on the genetic diversity of within a single beekeeping operation, apiary and honeybee colony, which to date remain poorly understood.

Methods

Outbreak description and microbiological examination

In Slovenia, an AFB zone is defined as a 3 km radius around the apiary with confirmed AFB and extended to a larger surveillance area with each newly detected AFB-positive apiary. In September–October 2019, AFB was confirmed by clinical and microbiological examination in 9/17 apiaries, which were maintained by six beekeepers and were located in a single AFB zone in northwestern Slovenia (Fig. 1, Table S1, available in the online version of this article). The remaining eight apiaries maintained by seven beekeepers showed no clinical symptoms of AFB (altered brood) and was detected only in honey samples. For microbiological examination, swabs of altered broods (symptomatic colonies) and/or minimum of 50 ml of honey (symptomatic and asymptomatic colonies) from individual colonies were both collected at clinical examination of apiaries by the veterinarian specializing in bee health.

Fig. 1.

Map of the investigated AFB outbreak. A total of 59 outbreak-related isolates were included in the study. The corresponding apiaries (1–17) are represented as circles; position corresponds to the geographical location of apiaries and colour to the identified outbreak clone. Apiaries 1–3 belonging to the same beekeeping operation (maintained by beekeeper 1) are written in red.

Isolation of was performed according to the recommendations of the World Organization for Animal Health [2]. Honey samples were pre-heated to 45–50 °C and ~50 g of liquefied honey was poured into sterile centrifuge bottles, supplemented with 150 ml of warm sterile water, mixed well to a homogeneous suspension and centrifuged for 30 min at 4000 . The supernatant was discarded, leaving 1–2 ml at the bottom of the centrifuge bottle including the pellet, which was then resuspended, transferred to a sterile sample tube and supplemented with 0.9% NaCl to a final volume of 10 ml. The suspension was incubated at 80 °C for 10 min to kill the vegetative cells. Brood swabs and 500 µl of the prepared honey suspensions were streaked onto brain heart infusion agar (Oxoid) plates supplemented with 5% sheep blood, 1 mg l–1 thiamine and 30 mg l–1 nalidixic acid, and incubated aerobically for 2–7 days at 37 °C. Suspect colonies with typical morphology were inoculated onto fresh plates to obtain in pure culture; after Gram staining and catalase testing, species identification was performed using MALDI-TOF MS (Bruker Daltonics). One colony per sample was selected, except for the three honey samples (3.1–3.3) from apiary 3 (Table S1) that served for determination of diversity within the colony and honey sample. The obtained isolates were stored at –70 °C and included in the strain collection of the Veterinary Faculty, National Veterinary Institute, Slovenia. The collection and the associated isolate metadata served for the selection of genetically related but epidemiologically unrelated (non-outbreak) isolates.

A total of 59 outbreak-associated isolates underwent ERIC-PCR and WGS typing (Table S1); of these, 37 originated from honey and 22 from a brood. The isolates originated from symptomatic and asymptomatic colonies from 17 apiaries maintained by 13 beekeepers. In total, 40/59 isolates originated from a single beekeeping operation comprising three geographically closely located apiaries (1–3) maintained by a single beekeeper (hereinafter referred to as ‘beekeeper 1’; Table S1), which enabled us to assess the genetic diversity of within a single beekeeping operation, apiary and colony (Fig. 2). For this, the sampling effort was increased in apiaries 1–3 maintained by the same beekeeper 1, totalling 40 isolates. Of these, 20 isolates originated from five colonies of apiaries 1 and 2 (five isolates from brood samples and five from honey samples) and 20 isolates from four colonies of apiary 3 (four isolates from brood samples and 16 from honey samples since for three honey samples, five isolates per sample were selected for further investigation) (Fig. 2, Table S1). Since 15 isolates originated from three honey samples, diversity within a honey sample and the corresponding colony (when interpreted in combination with the brood isolate) could be assessed.

Fig. 2.

Study overview of 59 outbreak-related isolates that underwent ERIC and WGS typing. For isolate metadata, see Table S1.

ERIC-PCR

ERIC typing was performed as previously described [5, 6], with minor modifications [14]. For ERIC pattern determination, amplicons were analysed by QIAxcel capillary electrophoresis (Qiagen) as described previously [14]. The prominent 2800 bp band, which is present in ERIC II isolates but absent from ERIC I isolates and has been described as important for the differentiation of ERIC types I and II [6, 7], was considered relevant in this study.

WGS and cg/wgMLST typing

Isolates were revived from frozen stocks (–70 °C) by streaking on blood agar plates, followed by 72 h of incubation at 37 °C. For ERIC-PCR and WGS, DNA was extracted using the DNeasy Blood and Tissue kit (Qiagen) according to the protocol for Gram-positive bacteria. The extracted DNA was quantified by Qubit 3.0 (Thermo Fisher Scientific) using the Qubit 1× dsDNA HS assay kit. DNA libraries were prepared using the NEBNext Ultra DNA Sample Prep Master Mix kit (New England Biolabs). Paired-end sequencing (2×150 bp) was performed on the NovaSeq 6000 system (Illumina) to a minimum coverage of 130×.

Raw reads were assembled using SPAdes version 3.7.1 implemented in BioNumerics version 7.6.3 (bioMérieux, Applied Maths NV). The wgMLST quality assessment window was used to evaluate the quality of the obtained reads and assemblies. The following thresholds were applied: (i) average read quality of Q>30 and expected coverage of >30 for the trimmed reads, and (ii) N 50 >21250 kb, number of contigs <380 and total assembly length of ~4.1 Mb (range=3.2–5.3 Mb) for the assemblies. All genomes passed the applied thresholds.

wgMLST analysis was performed in BioNumerics by applying the stable wgMLST scheme for consisting of 5745 loci [14]. Assembly-based and assembly-free allele calling were performed using default parameters. The UPGMA (unweighted pair group method with arithmetic mean) tree and the minimum spanning tree (MST) based on the allele profiles of 5738 wgMLST loci (excluding the partial sequences of seven MLST loci) were constructed in BioNumerics with the allele calls considered as categorical values. For easier comparison with previous cgMLST-based studies on AFB outbreaks [13, 15], a cgMLST analysis based on 2833 loci [14] was also performed.

Six genetically most closely related but epidemiologically unrelated (non-outbreak) isolates from the national WGS database originating from Slovenia (n=5) and Germany (n=1) were added to the analysis to assess the ability of WGS to differentiate between outbreak-related and non-related isolates, and to establish a threshold with the highest specificity (Table S1). The constructed UPGMA wgMLST tree was visualized and annotated using iTol version 4.4.2 [18].

In silico seven-gene MLST

MLST sequence types (STs) were extracted from the assembled genomes of 59 isolates using the MLST plugin implemented in BioNumerics; the PubMLST nomenclature [19] on seven previously described housekeeping genes [8] was used.

Visualization of outbreak data

The geographical location of each apiary was used to show the spatial distribution of the obtained isolates (from one to 20 isolates per apiary). Microreact [20] was used for visualization.

Results

For the investigated outbreak, both cgMLST and wgMLST showed the existence of two outbreak clusters, which could be delineated by applying a single-linkage threshold of 27 allele differences (AD) by cgMLST (Fig. S1) and 34 AD by wgMLST (Fig. 3). Because wgMLST showed a higher discriminatory power and better specificity to exclude the non-outbreak isolates as compared with cgMLST (see isolate PL189 in Figs S1 and 3), hereafter only the wgMLST results are discussed.

Fig. 3.

wgMLST minimum spanning tree of the 59 outbreak-related and six non-outbreak-related isolates. Isolates are coloured according to their outbreak association. Numbers on the connecting lines indicate allele differences (AD); branch lengths are scaled logarithmically. For clarity, some distances within outbreak cluster 1 are not shown, but were ≤11 AD. Isolates with ≤34 AD are shaded in grey.

Both identified outbreak clusters belonged to the ERIC II-ST11 type. Cluster 1 comprised 55 isolates from 15 different apiaries, including all 40 isolates originating from a single beekeeping operation (beekeeper 1; Table S1, Figs 3 and 4). No evident clustering was observed with respect to the beekeeper, sample type or apiary, further confirming that the cluster 1 isolates are clonal (Fig. 4). The median pairwise distance between isolates of cluster 1 was 11 AD (range=1–33) and between isolates of the same beekeeping operation (i.e. beekeeper 1) it was 10 AD (range=1–22) (Tables 1 and S2). The maximum pairwise AD between cluster 1 isolates was 33 (Table 1), and the isolates could be linked by a single-linkage distance of 16 AD on the wgMLST MST (Fig. 3). For cluster 1, the minimum distance between the outbreak-related and non-outbreak isolates was 39 AD. With respect to beekeeper 1 with three inspected AFB-positive apiaries and three colonies of apiary 3 inspected in detail, the median pairwise distances for the apiaries ranged from 8 to 11 (min.=1, max.=22) and for the colonies of the same apiary from 8 to 14 (min.=1, max.=20), regardless of whether brood isolates were included or not (Table 1).

Fig. 4.

UPGMA wgMLST tree of 59 outbreak-related and six non-outbreak isolates. Isolates are coloured according to their outbreak association, sample type and apiary (beekeeper 1). *Honey sample 3.1; **honey sample 3.2; ***honey sample 3.3 (all honey samples denoted by asterisks originated from apiary 3 maintained by beekeeper 1).

Table 1.

Median, minimum and maximum wgMLST allele differences (AD) between outbreak-associated isolates at different levels

Level	Median AD	Min. AD	Max. AD	No. of isolates
Outbreak cluster 1	11	1	33	55
Outbreak cluster 2	30	4	44	4
Beekeeping operation (beekeeper 1)	10	1	22	40
Apiary 1	8	1	14	10
Apiary 2	11	3	18	10
Apiary 3	10	1	22	20
Colony 3.1	9	3	13	6
Colony 3.2	8	1	16	6
Colony 3.3	14	1	20	6
Honey 3.1	8	3	13	5
Honey 3.2	10	4	16	5
Honey 3.3	14	1	20	5

Cluster 2 included four isolates from three different apiaries (Table S1, Fig. 3), which were separated from the isolates of cluster 1 by a minimum of 63 AD (data not shown). The median pairwise distance between cluster 2 isolates was 30 AD (range=4–44). The maximum pairwise AD between cluster 2 isolates was higher than for cluster 1 isolates, as well as the single-linkage distance (34 AD) (Fig. 3). For cluster 2, the minimum distance between the outbreak-related and non-outbreak isolates was 37 AD. Cluster 2 isolates formed a paraphyletic cluster with the non-outbreak isolate PL189 on the UPGMA wgMLST tree (Fig. 4), but clustered to the exclusion of PL189 on the wgMLST MST (Fig. 3). Because no available epidemiological evidence linked isolate PL189 to the studied outbreak and because it originated from a geographically distant (>80 km) apiary with respect to the outbreak-related isolates with no common foraging area between them, this isolate was regarded as a non-outbreak isolate.

An epidemiological link due to trade of honeybees and shared beekeeping activities (beekeeping equipment or beekeepers) was identified between apiaries 4 and 16, which were geographically distant (>3 km) and separated by various landscape features (i.e. river and landscape elevation), explaining the transmission of outbreak clone 2 between them (Fig. 1). We also identified an epidemiological link between several beekeepers who maintained different apiaries harbouring outbreak clone 1; namely, they were members of the same beekeepers’ association and exchanged beekeeping equipment. In addition, examples of inappropriate beekeeping practices (e.g. robbery of honeycombs) were reported in the investigated AFB zone, which contributed to the spread of AFB.

Isolates PL92 (belonging to cluster 2) and PL105 (belonging to cluster 1) both originated from apiary 4 maintained by beekeeper 2 (Fig. 1). Since the isolates of clusters 1 and 2 were spatiotemporally or epidemiologically related (Fig. 1), this suggests a concurrent AFB outbreak caused by two ERIC II-ST11 outbreak clones (biclonal outbreak).

Discussion

In this study, we investigated the genetic diversity of within a geographically restricted AFB outbreak and a single beekeeping operation. Using wgMLST, we identified two outbreak clusters, which could be delineated by applying the single-linkage threshold of 34 alleles, suggesting a concurrent outbreak with two clones. The isolates from a single beekeeping operation (≤22 pairwise AD) fell within outbreak cluster 1 (≤33 pairwise AD); diversity within the operation had already been reached when inspecting individual apiaries. A fairly high genetic diversity of was observed at all inspected levels (operation, apiary, colony and honey sample), which should be taken into account when interpreting the WGS results. The transmission events within the outbreak could be explained by the activities of honeybees (<3 km) and beekeepers.

Both outbreak clusters could be delineated by applying a single-linkage threshold of 34 and 27 AD on the wgMLST and cgMLST MSTs, respectively. These thresholds are similar to the thresholds of 35 and 24 AD, respectively, which we previously proposed [14]. By contrast, the 24-allele cgMLST threshold is markedly higher than the 10-allele threshold, which was based on the investigation of AFB outbreaks in Gotland (Sweden) using Ridom SeqSphere+ and different cgMLST schemes (ad hoc ERIC type-specific schemes and a stable 2419-loci scheme) [13, 15]. We re-analysed these outbreaks using BioNumerics and a 2833-loci cgMLST scheme (also used in this study) and found high agreement in terms of genetic relatedness. Thus, when applying the 10-allele threshold, all cgMLST-based approaches produce highly concordant results with regard to the observed genetic distances and cluster delineation. However, we argue that the cgMLST threshold is closer to 30 than to 10 AD since we set this threshold to 24 AD in our previous work [14] and to 27 AD in this study. When the 10-allele threshold was applied to the Gotland isolates [13], we observed the same clustering, but when this threshold was increased to 24 AD, the two closely related ST11 groups formed a single cluster (separated by ≤20 AD), which was supported by the underlying epidemiological data [13, 14]. The wgMLST analysis based on a 5738-loci scheme (also used in this study) showed that these two groups differed by ≤27 AD, further supporting their close relatedness [14]. Furthermore, when investigating the AFB outbreaks in Slovenia in 2017–2019, we showed that the outbreak-associated clonal representatives can differ in up to 24 cgMLST alleles [14].

The present study confirms that the threshold can exceed 10 AD by showing that clonal representatives can differ by up to 20 AD within a single honey sample or colony and 22 AD within a single apiary or beekeeping operation. Because an optimal threshold for cluster delineation should aim to encompass the highest possible number of clonal subpopulations within a single source, the present results suggest that such a threshold is probably closer to 34 and 27 AD, as determined here by wg/cgMLST MST analysis. These thresholds are also supported by the closely related non-outbreak isolates, which differed from the outbreak-related isolates in >34 wgMLST alleles and >23 cgMLST alleles in this study (Figs 3 and S1).

Of note, thresholds for cluster delineation should be used as flexible guidelines rather than absolute rules, and transmission links should be interpreted in the context of WGS tree topology and the underlying epidemiological data [10-12]. In the present study, the importance of epidemiological data and the use of flexible thresholds is exemplified by isolate PL189, which differed in a minimum of 24 cgMLST alleles and 37 wgMLST alleles from the outbreak-related isolates, but because no available epidemiological data linked the isolate to the investigated outbreak, this isolate was classified as of non-outbreak. This case also suggests that wgMLST can be used for additional confirmation of presumable outbreak clusters identified by cgMLST due to its higher discriminatory power, as previously described [14, 21, 22]. Different WGS analytical approaches, pipelines and pipeline settings can lead to different estimates of genetic relatedness and even different conclusions regarding outbreak cluster delineation when fixed thresholds for genetic relatedness are applied [10, 21]. In this study, a stable cg/wgMLST scheme with a standardized nomenclature was used to minimize the variation that occurs in the WGS data analysis step.

As already stated above, epidemiological data should always be considered when assessing the clonality of epidemiologically linked isolates (e.g. isolates originating from the same apiary or beekeeping operation). It is also important to distinguish between a single-linkage threshold based on MST and a maximum pairwise SNP/allele threshold based on pairwise distance matrix since the former tends to be lower than the latter, especially when a higher number of isolates is analysed. When the delineation of outbreak clusters is unclear and high discriminatory power is required, cgMLST can be further investigated by wgMLST and/or SNP approaches [14, 21, 22].

In this study, most transmission events could be explained by the activities of honeybees because apiaries were spatiotemporally linked (i.e. located in a single AFB zone). In addition, transmission of the outbreak clone between the three apiaries of the same beekeeping operation could also have occurred through the activities of beekeeper 1. An epidemiological link due to trade of honeybees and exchange of beekeeping equipment was identified for cluster 2 isolates, which supports the transmission between geographically distant apiaries.

In conclusion, this study reveals the genetic diversity of outbreak-causing clones at the outbreak, beekeeping operation, apiary, colony and honey sample levels. The genetic diversity of representatives within a single beekeeping operation, which reached 22 AD in this study, should be considered when assessing their clonality. The possibility of multiclonal outbreaks should not be neglected.

Click here for additional data file.

Click here for additional data file.

Click here for additional data file.