Mosquitoes provide a transmission route between possums and humans for Buruli ulcer in southeastern Australia

Ethics

Ethical approval for the use in this study of de-identified human Buruli ulcer case location, aggregated at the mesh-block level, was obtained from the Victorian government Department of Health Human Ethics Committee under HREC/54166/DHHS-2019-179235(v3), ‘Spatial risk map of Buruli ulcer infection in Victoria’.

Study site

Insects were collected from the Mornington Peninsula suburbs of Rye (population 8,416), Blairgowrie (population 2,313), Tootgarook (population 2,869) and Capel Sound (population 4,930)47 (Extended Data Fig. 1)48.

Arthropod collection and identification

Mosquito-trapping campaigns used Biogent Sentinal (BGS) traps (Biogent) that were baited with dry-ice pellets to provide a source of CO2 over an intended 12 hour period. BGS traps were set out at dusk, and collected at dawn, in shaded locations on the grassed edges of the roadways in the study area. GPS locations for all traps were recorded, with data collection managed using Atlas of Medical Entomology (v.3.4.4; Gaia Resources). Trapped mosquitoes were knocked down with CO2 by placing the catch bag in dry ice before being transported back to the laboratory and kept at −20 °C until processing. Mosquito species were morphologically identified using a stereo dissecting microscope (SMZ800N; Nikon) and reference to taxonomic keys49,50,51. Data collection was managed using Microsoft Excel (v.16.73).

To assess the potential presence of M. ulcerans beyond mosquitoes, other arthropods were collected using YST and SO. Two YST and SO were placed in residential properties where householders had previously noted insect activity to the researchers8. The SO were placed on the ground and had hay grass infusion (3 Jack Rabbit (clover or lucerne) pellets (Laucke Mills) in 500 ml water) added to them, whereas the YST (Bugs for Bugs) were placed on the ground with a 14 cm plant tag plastic T-support (Garden City Plastics). Within 3–4 days of being set, residents were asked to pack up the YST and SO by covering the sticky card with a plastic film and return the sealed traps to the laboratory, where they were stored at −20 °C. Non-mosquito arthropods were morphologically identified to family level and, if PCR positive for M. ulcerans, were DNA barcoded for species confirmation by targeting the COI gene52.

DNA extraction from mosquitoes

Mosquitoes were sorted by species per trap and by sex and then pooled in 2 ml O-ring tubes, with a maximum of 15 individuals in each pool. A subset of Ae. notoscriptus mosquitoes were also screened individually. Mosquitoes were homogenized with 10 mm × 1.0 mm zirconia/silica beads (BioSpec Products), with 597 µl of Buffer RLT and 2.8 µl of carrier RNA. Homogenization was performed using a TissueLyser II (Qiagen) at 30 oscillations per second for 100 s, repeated twice. Tubes were then centrifuged at 16,000g for 3 min. A 550 µl volume of supernatant was transferred into a 96-well deep-well plate, with extraction performed according to the protocol for the BioSprint 96 One-For-All Vet Kit (Qiagen). Every 11th of 12 wells in a 96-well plate was a blank DNA extraction control (7 in total) and a synthetic IS2404-positive control was spiked into 1 of these 7 wells to act as a positive extraction control. Extraction was performed on a KingFisher Flex Magnetic Particle Processor (Thermo Scientific).

DNA extraction from arthropods

Arthropods other than mosquitoes collected on YST or SO were removed from the sticky cards and placed in 1.5 ml microtubes (Eppendorf). Insects were separated by family and trap location, and were pooled with a maximum of 10 individuals from each family per 1.5 ml tube. Samples were extracted non-destructively to allow species confirmation if positive detections occurred. DNA was extracted using the ISOLATE II Genomic DNA Kit (Bioline). Briefly, 25 µl of Proteinase K and 180 µl of Lysis Buffer GL was added to each tube with samples incubated overnight at 56 °C. Following incubation, the insects were removed and stored to allow for further morphological identification if required, with the DNA extraction completed on the incubation solution as per manufacturer instructions.

Synthetic PCR positive control

A synthetic PCR positive control DNA molecule was designed to discriminate false positives due to contamination with positive control DNA versus the authentic IS2404 amplicon. The synthetic positive control was designed to have an amplicon size of 120 bp to easily differentiate it from a true IS2404 PCR positive (59 bp)15. The synthetic positive was added at the DNA extraction stage on all 96-well plates as a positive control for this step and for the subsequent qPCR. The additional DNA sequence used to construct the synthetic positive control was randomly selected from a DNA sequence unlikely to be in the laboratory, in this case, Irrawaddy dolphin (MK032252).

The synthetic positive control had the sequence 5′-TCCTAAAGCACCACGCAGCATCTATCGCGAGCTTAATCACCATGCCGCGTCCAACGCGATCCCCGCTCGGCAGGGATCCCTCTTCTCGCACCGGGCCACAATCCACTGGGGTCGCTATGA-3′ and was synthesized as a single-stranded DNA oligonucleotide (Sigma-Aldrich). The underlined regions indicate the forward primer, probe and reverse primer sequences, respectively. The synthesized IS2404 synthetic positive was resuspended in nuclease-free water and diluted to 0.001 pM, with 2 µl being used for extraction and positive controls. To confirm the presence of a true positive as opposed to contamination, 5 µl of the qPCR product was added to 1 µl of DNA Gel Loading Dye 6X (Thermo Scientific) and run on a 2% agarose gel (TopVision Agarose Tablets; Thermo Scientific), with 1% SYBR Safe DNA Gel Stain (Invitrogen). The size of any positive IS2404 detection was assessed against 2 µl of 100 bp DNA Ladder (Promega) and run at 50 V for 1.5 h before being visualized with an EZ Gel Documentation System (Bio-Rad). Before the screening of insects began, the synthetic positive control for IS2404 was successfully designed and tested. By running the amplified PCR products on an agarose gel with the synthetic positive control and a real positive control, visual differentiation could be determined between a synthetic positive occurring at 120 bp and a true positive at 59 bp (Supplementary Fig. 2).

Screening insects by qPCR for M. ulcerans

The qPCR screening was performed using three independent assays IS2404, IS2606 and KR15. All samples were first screened with the IS2404 qPCR; if a positive was detected, additional confirmation was attempted with IS2606 and KR qPCR assays. Reactions were performed using 7.5 µl TaqMan Fast Universal PCR Master Mix (2X), no AmpErase UNG (Applied Biosystems), 1 µl of the primer–probe mix, 2 µl DNA and 4.5 µl nuclease-free water. A final primer–probe concentration for the IS2404 assay was as follows: 250:650:450 nM for the forward primer, reverse primer, and probe; and 800:800:220 nM for the forward primer, reverse primer, and probe for the IS2606 and KR assay. A 2 µl volume of the synthetic positive control was added for the IS2404 reactions, whereas 2 µl of M. ulcerans DNA was used for IS2606 and KR. All reactions included 6 no-template extraction controls and were run in a 96-well-plate format. Cycling conditions were as follows: denaturation at 95 °C for 2 min, followed by 45 cycles at 95 °C for 10 s and 60 °C for 30 s, with qPCR performed on a QuantStudio 5 Real-Time PCR System (Applied Biosystems). Positives were classified as reactions that produced a cycle quantification value less than 40. Data were analysed using the QuantStudio Design and Analysis Software v.1.4.3 with the ΔRn threshold set at 0.04 for IS2404 and IS2606, and a ΔRn threshold of 0.1 for KR. The MLE per 1,000 mosquitoes tested (bias-corrected MLE for point estimation of infection rate and a skew-corrected score CI) was calculated from the pooled samples53. Fisher’s exact test for assessing the significance of differences in IS2404 PCR positivity between mosquito species was calculated in R 4.0.2 (ref. 54).

All qPCR screening was performed blind with mixed-species 96-well plates. The synthetic positive control was added at the DNA extraction phase to one well of each plate and to qPCR plates to check that both extraction and qPCR detections were successful. All no-template controls (extraction and qPCR stage) were checked to ensure they remained negative, and that synthetic positive controls were detected for both the DNA extraction and qPCR stage in each run. Positive samples were run on agarose gels to confirm they were true positives and not contamination from the synthetic positive control.

An IS2404 qPCR standard curve was prepared using 10-fold serial dilutions of M. ulcerans genomic DNA, with quadruplicate testing of each dilution. The DNA was extracted from M. ulcerans JKD8049 and quantified using fluorimetry (Qubit dsDNA HS; Thermo Fisher Scientific)20. A limit of detection was defined as the lowest dilution that returned a positive signal for all four replicates. GE were calculated based on the estimated mass of the M. ulcerans genome of 5.7 fg (ref. 20). IS2404 Ct values were converted to GE to estimate bacterial load within a sample by reference to a standard curve (r2 = 0.9956, y = [−3.829ln(x) + 37.17]Z, where r2 = correlation coefficient y = Ct and x = amount of DNA (in fg) and Z = the dilution factor; Supplementary Fig. 2). An IS2404 qPCR standard curve was fitted using nonlinear regression in GraphPad Prism (v.9.5.1) (Supplementary Fig. 3).

M. ulcerans genome sequencing

Whole-genome sequencing was performed directly on DNA extracted from selected PCR-positive mosquito samples and possum excreta specimens using a hybridization capture approach, based on 120 nt RNA baits spanning the 5.8 Mbp chromosome of the M. ulcerans JKD8049 reference genome (BioProject ID PRJNA771185) (SureSelect Target Enrichment System; Agilent; Supplementary dataset 4) and the Illumina Nextera Flex for Enrichment with RNA Probes protocol55. Resulting sequence reads were submitted to National Centre for Biotechnology Information (NCBI) GenBank and are available under BioProject PRJNA943595 (Supplementary Table 2).

M. ulcerans SNP calling, SNP imputation and phylogenetic analysis

To compare genomic variations between M. ulcerans clinical isolate genome sequences and sequence-capture enrichment datasets, we mapped the sequence reads and called nucleotide variations using Snippy (v.4.4.5) against a finished M. ulcerans reference chromosome, reconstructed from a Victorian clinical isolate (JKD8049; GenBank accession NZ_CP085200.1; https://github.com/tseemann/snippy). Although standard parameters, including a minimum coverage of 10×, were used for the clinical isolates and two possum sequence-capture datasets, the mosquito sequence-capture datasets had lower read coverage, necessitating the adjustment of parameters. Thus, the minimum coverage threshold was lowered to 1× to facilitate SNP calling for the mosquito sequence-capture datasets. The resulting SNPs were combined with 117 SNPs obtained from a reference set of 36 M. ulcerans genomes that represented the previously defined population structure of the pathogen in Victoria24. Due to the low read coverage however, the number of core variable nucleotide positions mapped among the 5 sequence-capture datasets was variable (range, 22–112 variable nucleotide positions). To enable inclusion of the sequence-capture datasets that had missing SNP sites, we used a multivariate imputation approach, using the IterativeImputer function from scikit-learn56. The combined alignment of 117 core genome SNPs from the sequence-capture datasets and clinical isolate genomes served as the foundation for inferring a maximum likelihood phylogeny. This phylogeny was established using the GTR model of nucleotide substitution and executed with FastTree (v.2.1.10)57. The incorporation of R packages phytools (v.1.0-1)58 and mapdata (v.2.3.1)59 allowed for the alignment of tree tips against a base map, facilitating the visualization of geographical origins of the samples. Further details, including the code used for missing SNP imputation and phylogeographical analysis, can be found in our GitHub repository60.

Ae. notoscriptus typing and species confirmation sequencing

Mosquito genotyping was performed by sequence comparisons of a partial fragment of the COI gene52. DNA was extracted using the above protocols. PCR was performed using 5 µl of 5× MyFi Reaction Buffer, 1 µl MyFi DNA polymerase, 5 µl DNA and primer concentrations27, with the reaction made up to 25 µl with nuclease-free water. Reaction conditions were as follows for COI: initial denaturation at 95 °C for 1 min, followed by 35 cycles at 95 °C for 20 s, 46 °C for 20 s and 72 °C for 60 s, before a final extension at 72 °C for 5 min. A 5 µl volume of the amplified PCR product was added to 1 µl of DNA Gel Loading Dye 6X (Thermo Scientific) and run on a 1% agarose gel (TopVision Agarose Tablets; Thermo Scientific), with 1% SYBR Safe DNA Gel Stain (Invitrogen). A 2 µl volume of 100 bp DNA Ladder (Promega) was added to confirm amplicon size and run at 100 V for 45 min. PCR products that produced bands of the correct size were purified using the ISOLATE II PCR and Gel Kit (Bioline), as per the manufacturer’s protocol, and submitted for sequencing using an Applied Biosystems 3730xl capillary analyser (Macrogen), with sequencing occurring on both strands. Sequences were analysed in Geneious Prime (v.2019.2.1) and trimmed to high-quality bases, aligned using ClustalW v.2.1 and trimmed to a consensus region, Ae. notoscriptus COI (874 bp) and for species identification COI (816–882 bp). Sequences were analysed using blastn against the NCBI database. COI sequences generated for species identification are available under accession numbers OQ600123–OQ6001234, and COI sequences for Ae. notoscriptus phylogenetics are under accession numbers OQ588831–OQ588867 (Supplementary Fig. 1).

Mosquito blood-meal analysis

Ninety blood-fed mosquitoes were identified as having an engorged abdomen and still having a red pigment (Sella score 2–3), indicating a fresh blood meal, and were dissected with a sterile scalpel blade. Blood from the dissected abdomen was absorbed onto a 3 mm × 20 mm piece of a Whatman FTA card (Merck) and placed in a 2 ml tube. DNA was extracted from the FTA card using an ISOLATE II Genomic DNA Kit (Bioline) with a pre-lysis in 180 µl of Lysis Buffer GL and 25 µl Proteinase K for 2 hours before completing the extraction as per the manufacturer’s protocol. Extracted DNA was amplified for cytb using primers previously described61, with MyTaq HS Red Mix (Bioline), and the thermocycling conditions used were as follows: 95 °C for 1 min, 30 cycles of 95 °C for 15 s, 50 °C for 20 s, 72 °C for 20 s and a final extension at 72 °C for 2 min. Negative extraction and negative PCR controls were included with each PCR reaction. A volume of 5 µl of the amplified products was run on a 1% agarose gel containing 0.1% SYBR Safe DNA gel stain (Invitrogen) and visualized with a G:BOX Syngene blue-light visualization instrument. If a band was visualized at ∼480 bp, the PCR product was purified using an ISOLATE II PCR and Gel Kit (Bioline).

PCR products were quantified using a Qubit (Invitrogen) fluorometer with an HS dsDNA kit. Sequencing libraries were prepared using 10 ng of purified PCR product with a NEXTFLEX Rapid XP DNA-Seq Kit (PerkinElmer) barcoded using NEXTFLEX UDI Barcodes (PerkinElmer). As a result, 70 blood-meal libraries were sequenced, along with 6 PCR negative controls and 3 extraction negatives, on a NovaSeq 6000 (Illumina), with 2 Gb requested per sample.

Sequence data were analysed by identifying poor-quality reads using Rcorrector62 and removed with TrimGalore v.0.6.5. De novo assembly was performed on the remaining sequence reads using Trinity v.2.8.663. The resulting contigs were filtered to be between 400 bp and 480 bp and analysed using BLASTN against the nucleotide database publicly available on the NCBI website. The resulting hits were filtered to exclude those sequences that had 64. Sequence reads were submitted to GenBank (BioProject ID PRJNA943595).

Mosquito phylogenetic analysis

Mosquito genotyping was performed with COI because this genetic marker has previously been used to identify potential cryptic species within Ae. notoscriptus and can provide better resolution than other markers such as ND5, CAD or EPIC (exon-primed intron crossing)27. Phylogenetic analysis was performed on trimmed consensus regions of Ae. notoscriptus COI. The substitution model was selected using jModelTest2 v.2.1.10, with the topology taking the best of nearest-neighbour interchange, subtree pruning and regrafting65. The most appropriate substitution model was selected based on the Akaike information criterion. Maximum likelihood trees were constructed in PhyML v.3.3.2 with 1,000 bootstrap replicates; the gamma distribution parameter was used to estimate rate variation across sites66. The Hasegawa–Kishino–Yano (HKY) substitution model was selected for the COI tree.

Geographical data acquisition and spatial cluster analysis

The population map was created using QGIS geographical information system software (v.3.16.7)67, using a 1 km2 population grid68 with 2011 Victorian mesh-block data. Since 2004, Buruli ulcer has been a notifiable condition in Victoria, requiring health department reporting by doctors and laboratories. De-identified case notification data of patients with Buruli ulcer who had laboratory-confirmed M. ulcerans infection and who lived on the Mornington Peninsula during the years 2019–2020 were provided by the Victorian Department of Health. The cases were defined as patients with a clinical lesion that was diagnosed using IS2404 qPCR and culture15. To conduct high-resolution spatial analyses, the data were aggregated at the mesh-block level, the smallest geographical census units which typically contain 30–60 dwellings. The 2011 Victorian mesh-block boundaries and the Victorian mesh-block census population counts datasets were obtained from the Australian Bureau of Statistics website69. The datasets were joined using the unique mesh-block IDs using QGIS (v.3.16.7)67. The latitude and longitude (projected in GDA94) were derived from the centroids of the mesh-block polygon. The dataset was then downsampled to include only the Mornington Peninsula study area, specifically the Point Nepean and Rosebud–McCrae Australian Bureau of Statistics level 2 Statistical Area.

SaTScan v.10.1.0 (ref. 70) was used to identify spatial clusters among trapped mosquitoes positive for M. ulcerans, possum excreta positive for M. ulcerans and human Buruli ulcer cases. The software searches for instances where the observed number of spatial incidences exceeds the expected number within a circular window of varying size across a defined study area. A log-likelihood ratio statistic is calculated for each window by comparing the number of observed and expected cases inside and outside the circle against the assumption of randomly distributed cases. In addition to the most likely cluster, there are usually secondary clusters with almost as high likelihood that substantially overlap with the primary clusters. These secondary clusters can be indicative of subclusters within the primary cluster or potentially distinct clusters that are spatially adjacent to the primary cluster. The Mornington Peninsula surveillance data used in these analyses consisted of trapped mosquitoes (177 traps screened for IS2404 collected 12 November 2019 to 20 March 2020), M. ulcerans detected in possum excreta collected during the summer (December to February) of 2019 using data from a previous study26 and notified human Buruli ulcer cases from the study area in the years 2019–2020. The use of possum excreta collected during the summer was appropriate as Buruli ulcer transmission is most likely to occur during that time of year71,72. For each of the three data sources (trapped mosquitoes, possum excreta and human cases), the null hypothesis assumes that M. ulcerans detections or Buruli ulcer cases are uniformly distributed across the study area, where the alternative hypothesis suggests that there may be certain locations with higher rates than expected if the risk was evenly distributed. Primary and secondary clusters were accepted only if the secondary clusters did not overlap with previously reported clusters with a higher likelihood. Given that the trapped mosquito and possum excreta IS2404 PCR results were binary (positive or negative), a Bernoulli model was used to scan for spatial clusters, with the maximum cluster size set to 50% of the population size. The human Buruli ulcer case data aggregated at the mesh-block level varied in number, with some mesh blocks having zero cases and others having one or more. We applied the Poisson probability model to the notified Buruli ulcer case counts, using a background population at risk that was derived from the 2011 population census. The maximum cluster size was limited to 14,481 individuals, 10% of the total population at risk. To determine the likelihood of a triple-cluster overlap between the three SaTScan analyses occurring by chance, we conducted a permutation test73. In each of 10,000 iterations, the geographical coordinates for each variable were randomly shuffled. The number of SaTScan clusters with triple overlap was determined using the sf package74 in the R statistical programming language54.

Reporting summary

Further information on research design is available in the Nature Portfolio Reporting Summary linked to this article.