3 Microbiome of ticks and wildlife

Preface

Attribution Statement

The following chapter has been drafted in accordance with the journal Microbial Genomics.

The following manuscript has been published: Egan, S., Taylor, C., Banks, P., Northover, A., Alhstrom, L., Ryan, U., Irwin, P., and Oskam C. 2021. The bacterial biome of ticks and their wildlife hosts at the urban-wildland interface. Microbial Genomics, 7(12): 000730. DOI: 10.1099/mgen.0.000730

The following authors contributed to this manuscript as outlined below⁴.

Authorship order	Contribution (%)	Concept Development	Data Collection	Data Analysis	Drafting of manuscript
Siobhon L. Egan	70.0	X	X	X	X
Casey L. Taylor	6.0		X
Peter B. Banks	2.5	X
Amy S. Northover	2.5		X
Liisa A. Ahlstrom	2.5		X
Una M. Ryan	5.0	X
Peter J. Irwin	5.0	X
Charlotte L. Oskam	6.5	X

By signing this document, the Candidate and Principal Supervisor acknowledge that the information provided is accurate and has been agreed to by all other authors.

__________________ __________________
Candidate Principal Supervisor

Chapter linking statement: This chapter looks at the bacterial microbiome of wildlife samples. Using amplicon metabarcoding targeting the 16S rRNA gene to characterise the suite of bacteria in blood, tick and tissue samples. It provides an unbiased approach to investigate (i) what microbes are present in these wildlife samples and (ii) uncover overlap between sample types to provide insight into potential bacteria that maybe transmitted into wildlife which act as reservoir hosts. To obtain more detailed genetic information on ‘taxa of interest’ further sequencing was performed. This chapter provides important insights into tick (and potentially other vector) associated microbes in Australia.

Acknowledgment statement: This work was supported by resources provided by the Pawsey Supercomputing Centre with funding from the Australian Government and the Government of Western Australia. We thank Dr David Chandler, Dr Christopher Noune and Dr Matthew Stevens at the Australian Genomics Research Facility for their assistance. We thank the Department of Biodiversity, Conservation and Attractions staff for their invaluable assistance in Western Australian fieldwork, in particular Rebecca Kay, Douglas Giles and Hannah Kilian. We also thank Jade Kelly and volunteers for assistance in fieldwork and sample collection.

Funding statement: This study was part-funded by the Australian Research Council (LP160100200), Bayer HealthCare (Germany) and Bayer Australia. S.L.E. was supported by an Australian Government Research Training Program (RTP) Scholarship. C.L.T. was supported by a scholarship from the Northern Beaches Council. This project was also part supported by The Holsworth Wildlife Research Endowment from The Ecological Society of Australia (awarded to S.L.E.) and the Paddy Pallin Science Grant from The Royal Zoological Society (awarded to C.L.T.).

Data availability: Raw Illumina MiSeq data has been deposited in the European Nucleotide Archive under project number PRJEB46006 (16S rRNA bacteria) and PRJEB46056 (12S rRNA ticks). Nucleotide data produced by Sanger sequencing has been submitted to NCBI GenBank nucleotide database under accession numbers MW633074, MW633076–MW633082, MW633160–MW633167. A subset of representative zOTUs generated from tick 12S rRNA metabarcoding have been submitted to GenBank under accession numbers MW665133–MW665150. Additional information is available within supplementary files including links to electronic files which have been deposited on FigShare repository. Code used for analysis and RData files used for bioinformatics and statistical analysis are available on GitHub repositories siobhon-egan/wildlife-bacteria and siobhon-egan/wildlife-ticks.

Animal ethics: Murdoch University Animal Ethics Committee permit number R3026/18. Department of Biodiversity Conservation and Attractions permit numbers 2018/54B, 2018/55B, 2018/57D and 2017/30. University of Sydney Animal Ethics Committee Permit number 2018/1429.

Author contributions:: Conceptualization: S.L.E., U.M.R., P.J.I., C.L.O. Data curation: S.L.E., C.L.T. Formal Analysis: S.L.E. Funding acquisition: S.L.E., C.L.T., P.B.B., L.A.A., U.M.R., P.J.I., C.L.O. Investigation: S.L.E., C.L.T. Methodology: S.L.E., C.L.T., A.S.N. Project administration: P.J.I., C.L.O. Resources: A.S.N., P.B.B., P.J.I., C.L.O. Software: S.L.E. Supervision: P.B.B., U.M.R., P.J.I., C.L.O. Visualization: S.L.E. Writing – original draft: S.L.E. Writing – review & editing: S.L.E., C.L.T., P.B.B., A.S.N., L.A.A., U.M.R., P.J.I., C.L.O.

Keywords: Tick, microbiome, wildlife, sylvatic cycle, tick-borne diseases,

3.1 Abstract

Advances in sequencing technologies have revealed the complex and diverse microbial communities present in ticks (Ixodida). As obligate blood feeding arthropods, ticks are responsible for a number of infectious diseases that can affect humans, livestock, domestic animals and wildlife. While cases of human tick-borne diseases continue to increase in the northern hemisphere, there has been relatively little recognition of zoonotic tick-borne pathogens in Australia. Over the past five years, studies using high-throughput sequencing technologies have shown that Australian ticks harbour unique and diverse bacterial communities. In the present study, free-ranging wildlife (n = 203), representing ten mammal species, were sampled from urban and peri-urban areas in New South Wales (NSW), Queensland (QLD), and Western Australia (WA). Bacterial metabarcoding targeting the 16S rRNA locus was used to characterise the microbiomes of three sample types collected from wildlife; blood, ticks and tissue samples. Further sequence information was obtained for selected taxa of interest. Six tick species were identified from wildlife: Amblyomma triguttatum, Ixodes antechini, Ixodes australiensis, Ixodes holocyclus, Ixodes tasmani and Ixodes trichosuri. Bacterial 16S rRNA metabarcoding was performed on 536 samples and 65 controls, generating over 100 million sequences. Alpha diversity was significantly different between the three sample types, with tissue samples displaying the highest alpha diversity (P < 0.001). Proteobacteria was the most abundant taxon identified across all sample types (37.3%). Beta diversity analysis and ordination revealed little overlap between the three sample types (P < 0.001). Taxa of interest included Anaplasmataceae, Bartonella, Borrelia, Coxiellaceae, Francisella, Midichloria, Mycoplasma and Rickettsia. Anaplasmataceae bacteria were detected in 17.7% (95/536) of samples and included Anaplasma, Ehrlichia and Neoehrlichia species. In samples from NSW ‘Ca. Neoehrlichia australis’, ‘Ca.. Neoehrlichia arcana’, Neoehrlichia sp. and Ehrlichia sp were identified. A putative novel Ehrlichia sp. was identified from WA and Anaplasma platys was identified from QLD. Nine rodent tissue samples were positive for a novel Borrelia sp. which formed a phylogenetically distinct clade separate from the Lyme Borrelia and relapsing fever groups. This novel clade included recently identified rodent-associated Borrelia genotypes, which were described from Spain and North America. Bartonella was identified in 12.9% (69/536) of samples. Over half of these positive samples were obtained from black rats (Rattus rattus), and the dominant bacterial species identified were Bartonella coopersplainsensis and Bartonella queenslandensis. The results from the present study show the value of using unbiased high throughput sequencing applied to samples collected from wildlife. In addition to understanding the sylvatic cycle of known vector associated pathogens, surveillance work is important to ensure preparedness for potential zoonotic spillover events.

3.2 Introduction

Ticks carry a diverse range of infectious microbes such as viruses, piroplasms, spirochaetes and Rickettsiales (Pfaffle et al. 2013). Globally recognised infectious human tick-borne diseases include Lyme borreliosis, ehrlichiosis, babesiosis, tick-borne encephalitis and Powassan viral disease. Additionally, tick bite can cause reactions such as paralysis and anaphylaxis (Van Nunen et al. 2009). The incidence of tick-borne diseases in many parts of the world are rapidly increasing in both prevalence and geographic distribution (Madison-Antenucci et al. 2020). The sylvatic cycles of tick-borne pathogens in the northern hemisphere are generally well understood with respect to competent tick vectors and the reservoir hosts (e.g. Babesia microti (Westblade et al. 2017) and Borrelia burgdorferi sensu lato (Radolf et al. 2012)). The value of wildlife health surveillance as a tool for the detection of emerging zoonotic infectious disease has been demonstrated in outbreaks historically, such as malaria (Rondon et al. 2019), and more recently with SARS-CoV-2 (Plowright et al. 2021).

In Australia, three human tick-borne diseases are recognised; Queensland tick typhus (Rickettsia australis), spotted fever (Rickettsia honei subspp.) and Q fever (Coxiella burnetii) (S. R. Graves and Stenos 2017). Over the past decade there have been increasing concerns about unidentified tick-borne illness (syn. Lyme-like disease) in Australia (Chalada, Stenos, and Bradbury 2016), recently designated ‘Debilitating Symptom Complexes Attributed To Ticks’ (DSCATT) by the Australian Government Department of Health (Radcliffe et al. 2016).

Presently, 74 tick species have been described in Australia, including five introduced species (D. Barker 2019; Evans et al. 2019), of which 22 species have been identified biting humans (Kwak 2018). However, data show that Ixodes holocyclus and Amblyomma triguttatum account for >90% of human tick bites, with the remainder of reports largely attributed to Ixodes tasmani, Haemaphysalis bancrofti, Haemaphysalis longicornis and Bothriocroton hydrosauri (Kwak 2018). Small to medium-sized mammals are important in maintaining these tick populations, in particular during the immature life stages (Murdoch and Spratt 2005; Lydecker, Hochuli, and Banks 2019; Taylor et al. 2020). While data on the incidence of human tick bites in Australia is limited, evidence suggests that tick bites are common in peri-urban areas, including residential yards (C Taylor, pers comms).

Over the past decade, due to new technologies there is mounting evidence that Australian ticks harbor a unique array of microbes (Gofton, Oskam, et al. 2015; Siobhon L. Egan et al. 2020). In addition to these molecular studies, more targeted screening of potential reservoirs such as dogs (Irwin et al. 2017) and sea bird ticks (Moon et al. 2018) have added to the growing body of evidence suggesting that human tick-borne pathogens described in the northern hemisphere are not endemic in Australia. Some of the microbes that have recently been identified from Australian ticks however include new species of Borrelia (S. M. Loh et al. 2017; Panetta et al. 2017), Ehrlichia (Gofton et al. 2017), Neoehrlichia (Gofton et al. 2016), Babesia, and Theileria (T. L. Greay, Zahedi, et al. 2018; S. M. Loh, Egan, et al. 2018) and several viruses (O’Brien et al. 2018; Harvey et al. 2019). Little is known about the sylvatic life cycle of these microbes or their wildlife species reservoir hosts. Even for the few tick-borne pathogens recognised in Australia (e.g. Rickettsia spp.), much of this information remains elusive.

For a tick-borne disease to become established it requires a temporal and spatial overlap between the microbe (pathogens), competent vector tick(s), and one or more vertebrate reservoir hosts. In the northern hemisphere studies have shown that there are important differences in the role of host species in maintaining the sylvatic cycle of tick-borne pathogens (Levi et al. 2016; Huang et al. 2019; Stewart Merrill and Johnson 2020). Broadly, reservoir hosts can be classified as amplification and dilution hosts. Amplification hosts, in most cases, are small mammals that host immature ticks (larvae or nymph stages) and increase the environmental density of infected nymphs (Levi et al. 2016). In contrast, dilution hosts are inefficient reservoirs of the microbial agent(s) (Takumi, Sprong, and Hofmeester 2019; Ginsberg et al. 2021). An additional factor is the role of hosts throughout the tick life cycle. Larger animals are considered important in maintaining adult stages of the tick; white-tailed deer (Odocoileus virginianus) and Ixodes scapularis (Ostfeld et al. 2018) for example. In the Australian context it is generally assumed that small mammals such as bandicoots, possums and rodents are likely reservoirs for tick-associated microbes. In turn, larger mammals such as kangaroos and wallabies may then play an important role in maintaining tick populations, filling a similar ecological niche as deer in the northern hemisphere. While studies have used molecular tools to investigate microbes present within ticks in Australia (Gofton et al. 2016; S. M. Loh, Egan, et al. 2018; Moon et al. 2018), very little research has been conducted at the vertebrate host level to further elucidate the sylvatic cycle of tick transmitted bacteria. Bacterial amplicon sequencing is valuable in the setting of wildlife epidemiology surveys where investigations are not confined to existing presumptions of microbes (Galan et al. 2016), however such studies of the ticks alone, without concurrent analysis of vertebrate host tissues, cannot distinguish infectious components of the tick’s total microbiome from those that are non-transmissible.
Therefore, the search for tick-associated microbes in vertebrate host(s) is the next logical piece of information required to untangle the complex dynamics of tick-borne pathogens. In addition to novel pathogen discovery, high throughput sequencing technologies applied concomitantly to tick and vertebrate host tissues can provide insights into the dynamics of microbial epidemiology and interactions within landscapes (Titcomb, Jerde, and Young 2019). Identification of microbes may inform our understanding of the tick-borne disease risk and potential differences among temporal and spatial gradients.

The present study used metabarcoding to explore the bacterial communities present in Australian wildlife and their ticks, with a focus on bacteria related to known (vector-borne) pathogens, to provide insights into the sylvatic life cycle of tick associated organisms. This research was conducted in an ecological zone referred to as the urban-wildland interface, where there is a transition between wilderness and land developed by human activity – an area where a built environment meets or intermingles with a natural environment.

3.3 Methods

3.3.1 Study sites and sample collection

Small mammal trapping

Small mammal trapping was performed at various sites in Perth, Western Australia and Sydney, New South Wales (Figure 3.1). Each site was targeted in urban and peri-urban areas to reflect the urban-wildland interface. Elliot and cage traps were set and baited with universal bait (rolled oats, peanut butter and sardines). Traps were set at dusk and cleared at sunrise over three or four consecutive nights. To ensure a comprehensive assessment and sampling, target mammals were briefly anaesthetised in the field using isoflurane (I.S.O. 1mL/mL) vaporised in medical oxygen. Animals were weighed and examined systematically for ectoparasites. Up to 1 mL of blood was collected from either the caudal (tail) vein, femoral vein or ear capillary and stored in Mini Collect EDTA tubes. If practical, a 2 mm punch biopsy was taken at the tick bite site; where possible, a biopsy was taken from the ear and stored in RNAlater or 80% ethanol. Animals were systematically examined for ectoparasites which were removed and stored in 80% ethanol. Animals were recovered from anaesthesia by providing medical oxygen, and once fully alert were released at trap point. An individual mark was applied to identify animals using either a microchip administered subcutaneously or a unique patch of hair was removed. The number of animals and samples analysed, and locations from the present study are available in Table 3.1 and Figure 3.1 respectively.

Opportunistic collection

In a small number of cases, animal carcasses were also obtained opportunistically. These specimens were collected through incidental findings during fieldwork or from situations where animals were humanely euthanised in accordance with animal ethics and state and federal guidelines.

Tick identification

Ticks were collected and stored in tubes containing 70–90% ethanol. Samples were visualised using an Olympus SZ61 stereomicroscope (Olympus, Centre Valley, PA, United States) with an external Schott KL 1500 LED (Schott AG Mainz, Germany) light source. Photographs of tick specimens were taken using an Olympus SC30 digital camera (Olympus, Centre Valley, PA, United States) and analysis getIT software (Olympus, Centre Valley, PA, United States). Instar, sex and species was identified using a combination of available morphological keys and species descriptions (F. S. H. Roberts 1970; J. Jackson et al. 2002; Laan, Handasyde, and Beveridge 2011a; Stephen C. Barker and Walker 2014; Kwak et al. 2017).

Table 3.1: Animal sampling numbers for bacterial microbiome sequencing. Number of individual animals () by host species with number of genomic DNA samples analysed in parentheses.
Common name	Scientific name	n (individuals)	Blood	Tick	Tissue
Black rat	Rattus rattus	88	68 (84)	54 (79)	71 (80)
Brown antechinus	Antechinus stuartii	5	0 (0)	4 (6)	2 (2)
Brush tail possum	Trichosurus vulpecula	27	18 (21)	15 (33)	16 (20)
Bush rat	Rattus fuscipes	3	2 (2)	3 (4)	3 (3)
Chuditch	Dasyurus geoffroii	22	22 (24)	1 (1)	9 (9)
Deer	Rusa timorensis	3	0 (0)	0 (0)	3 (6)
Long-nosed bandicoot	Perameles nasuta	44	20 (25)	31 (60)	40 (44)
Quenda	Isoodon fusciventer	2	1 (1)	2 (2)	0 (0)
Rabbit	Oryctolagus cuniculus	7	0 (0)	7 (20)	7 (7)
Swamp rat	Rattus lutreolus	2	1 (1)	0 (0)	1 (1)
Total		203	132 (159)	117 (205)	152 (172)

Map of study sites for collection of wildlife samples used in bacterial profiling. Sampling sites denoted by blue triangles and capital cities by pink circles (for geographical reference). Insert maps of sites in (1) Perth, Western Australia and (2) Sydney, New South Wales.

Figure 3.1: Map of study sites for collection of wildlife samples used in bacterial profiling. Sampling sites denoted by blue triangles and capital cities by pink circles (for geographical reference). Insert maps of sites in (1) Perth, Western Australia and (2) Sydney, New South Wales.

3.3.2 DNA Extractions

Total genomic DNA (gDNA) was extracted from 200 $\mu$L of blood using a MasterPure DNA purification kit (EpicentreBiotechnologies, Madison, Wisconsin, U.S.A) following the manufacturer’s recommendations. Where 200 $\mu$L of blood was not available, sterile DNA free phosphate-buffered saline (PBS) was used to make samples up to 200 $\mu$L. Genomic DNA was eluted in 30–40 $\mu$L of TE buffer and stored at -20$^\circ$C until further processing.

Tissue samples (skin and spleen) were first rinsed in sterile, DNA-free PBS and cut into small pieces (<1-2 mm) using a sterile scalpel blade. Samples were homogenised in 180 $\mu$L of buffer ATL and 20 $\mu$L of proteinase K was added and incubated at 56$^\circ$C for ~ 16 hours. gDNA was extracted using the QIAamp DNA Mini Kit (QIAGEN, Germany) following the manufacturer’s protocols with the exception that final elution volume was decreased to 40–50 $\mu$L to increase gDNA yield.

Once ticks were identified, they were surface-sterilised by washing in 10 % hypochlorite solution, rinsed in 70 % ethanol and DNA-free PBS, and then air-dried. Genomic DNA was extracted using the DNeasy Blood and Tissue kit (QIAGEN, Germany) for adults, or the QIAamp DNA Mini Kit (QIAGEN, Germany) for nymphs and larvae. Due to the large number of immature tick stages collected from some animals and the expected low DNA yield, up to 5 specimens were pooled for extraction based on host, instar, engorgement status and species identification as determined by morphological methods. Ticks were placed in a 2 mL safe lock Eppendorf tube with a 5 mm steel bead, frozen in liquid nitrogen for 1 min and homogenised by shaking at 40 Hz in a Tissue Lyser LT (QIAGEN, Germany). Final elution of DNA was adjusted according to tick size, between 30–150 $\mu$L of AE buffer was added to the silicon membrane.

Extraction controls (EXBs) consisting of 200 $\mu$L sterile DNA-free PBS were included randomly in each extraction batch (total EXBs = 34). Blood and tissue samples were extracted in a separate laboratory, away from where ticks were processed, to avoid cross-contamination between sample types. Sterile procedures were followed throughout the laboratory process.

3.3.3 Metabarcoding sequencing

A high-throughput metabarcoding approach was used to sequence samples using the Illumina MiSeq platform. Libraries were built following the 16S Metagenomic Sequencing Library Preparation (Illumina Part # 15044223 Rev. B), with amplicon PCR primers containing Illumina MiSeq adapter sequences (Table 3.2). A total of 536 samples from 203 individuals (159 blood, 205 tick pools, and 172 tissue samples) underwent bacterial profiling.

Bacterial 16S rRNA libraries were built targeting the 16S hypervariable region 1-2 using primers 27F-Y and 338R (Gofton, Oskam, et al. 2015). Reactions were carried out in 25 $\mu$L volumes each containing: 1X buffer (KAPA Biosystems, South Africa), 1.5 mM MgCl₂, 0.4 mg/mL of bovine serum albumin (Fisher Biotech), 0.4 $\mu$M of each primer, 0.25 mM of each dNTP, 0.5 U of Taq (KAPA Biosystems, South Africa) and 2 $\mu$L of gDNA. Thermal cycling conditions were as follows; 95$^\circ$C for 5 mins, followed by 35 cycles of 95$^\circ$C for 30 secs, 55$^\circ$C for 30 secs, 72$^\circ$C for 45 secs; and a final extension of 72$^\circ$C for 5 mins. For a subset of samples where microbes of interest were identified after initial screening at the v1-2 hypervariable region, additional bacterial 16S rRNA libraries targeting the v3-4 hypervariable region were also prepared. Primers 338F and 806R were used to target an ~450 bp product. Reactions and thermal cycling conditions were carried out as per the 27F-Y / 338R assay, except for an increased concentration of MgCl₂ to 2.0 mM.

To confirm morphological identification tick gDNA underwent amplicon metabarcoding approach targeting the mitochondrial 12S rRNA gene. Pan-Ixodida primers T1B/T2A (Beati and Keirans 2001) were used to amplify an ~370 bp product of the 12S rRNA Reactions were carried out as per the bacteria 27F-Y / 338R assay, except only 1 $\mu$M of tick gDNA template was added. Thermal cycling conditions were as follows; 95$^\circ$C for 5 mins, followed by 5 cycles of 95$^\circ$C for 15 secs, 51$^\circ$C for 30 secs, 68$^\circ$C for 30 secs; 25 cycles of 95$^\circ$C for 30 secs, 53$^\circ$C for 30 secs, 72$^\circ$C for 1 min; and a final extension of 72$^\circ$C for 5 mins.

Amplicon PCR products were then indexed using the Nextera XT DNA library preparation kit in 25 $\mu$L volumes following the manufacturer’s recommendations. All PCRs included no-template controls (NTC; total = 31) during each reaction set up and PCRs were performed under strict laboratory conditions. Amplicons were then dual-indexed using the Nextera XT index kit. Reactions were performed in 25 $\mu$L volumes following manufacturers recommendations. Libraries were purified with Axygen PCR clean up beads and quantified using Qubit High Sensitivity dsDNA assay kit (Thermo Fisher Scientific, Waltham, MA, USA) and pooled in equimolar amounts. Libraries were shipped to the Australian Genomic Research Facility (Melbourne, Australia) for final QC and sequenced on an Illumina MiSeq using v2 chemistry (2 x 250 paired-end).

Table 3.2: List of sequencing used for metabarcoding and target PCRs to screen blood, tissue and tick samples for bacteria. References 1. Gofton et al. (2015), 2. Turner et al. (1999), 3. Lopez et al. (2003), 4. Caporaso et al. (2011), 5. Beard et al. (2011), 6. Anderson et al. (1991), 7. Paddock et al. (1997), 8. Kawahara et al. (2004), 9. Beati and Keirans 2001 (2001).
Primer	Sequence (5’-3’)	Ref
Bacteria 16S rRNA
27F-Y	AGAGTTTGATCCTGGCTYAG	1
338R	TGCTGCCTCCCGTAGGAGT	2
338F	ACTCCTACGGGAGGCAGCAG	3
806R	GGACTACHVGGGTWTCTAAT	4
Bartonella 16S rRNA/ ITS
438s	GGTTTTCCGGTTTATCCCGGAGGGC	5
1100as	GAACCGACGACCCCCTGCTTGCAAAGC	5
Anaplasmataceae 16S rRNA
EC9	TACCTTGTTACGACTT	6
EC12	TGATCCTGGCTCAGAACGAACG	7
A171a	GCGGCAAGCCTCCCACAT	8
IS58-1345r	CACCAGCTTCGAGTTAAACC	8

3.3.4 Bioinformatics

Metabarcoding 16S rRNA sequence data was analysed using Quantitative Insights into Microbial Ecology (QIIME 2 2020.11) (Bolyen et al. 2019). Briefly, raw sequences were demultiplexed and quality filtered using the q2-demux plugin, followed by denoising using DADA2 (via q2-dada2) (Callahan et al. 2016) resulting in amplicon sequence variations (ASVs) (i.e. 100% identical operational taxonomic units or OTUs (Callahan, McMurdie, and Holmes 2017)). Bacterial data taxonomy was assigned to ASVs using the q2‐feature‐classifier (Bokulich et al. 2018) classify‐sklearn naive Bayes taxonomy classifier against the SILVA database (Quast et al. 2013) (release 132). Taxonomic assignments were confirmed using BLAST analysis (BLASTN 2.11.0+ (Zhang et al. 2000; Morgulis et al. 2008)) against the NCBI nucleotide collection (nt) database (accessed January 2021). Taxonomic lineage was then retrieved from NCBI taxonomy database using TaxonKit (Weiss et al. 2021) and adjusted to the lowest common ancestor based on percentage identity and e-value score.

Tick 12S rRNA metabarcoding sequences generated were analysed USEARCH v11 (R. C. Edgar 2010) pipeline. Briefly paired end reads were merged and sequences matching forward and reverse primers were retrieved (max number of mismatches = 2). Sequences were then quality filtered and singletons were removed. The unoise3 (Robert C. Edgar 2016) algorithm was used to perform denoising (error-correction) and generate zero-radius taxonomic units (zOTUs) equivalent to amplicon sequence variants (ASVs). Taxonomy was assigned using BLAST analysis and lineage retrieved as outlined above for 16S rRNA metabarcoding.

Data visualization and statistical analysis was performed in RStudio (v1.4) with R version 4.0.2. The main R packages used were ampvis2 v2.6.7, phyloseq v1.34.0, microbiomeutilities v1.00.12 and vegan v2.5-7. The R package decontam v1.10.0 (Davis et al. 2018) was used to inspect data for cross-contamination and cross-talk between samples and to establish read cut-off thresholds (prevalence threshold = 0.05). Identified contaminant ASVs (n = 104), controls samples (n = 65) and sequence values < 100 were removed. Rarefaction curves were generated using a step size of 100. Alpha diversity was measured using four indexes (number of observed ASVs, Chao1, Shannon and invSimpson) and statistical analysis was calculated using Wilcoxon pairwise test (non-parametric) between sample types (blood, tick and tissue). Heatmaps were generated using ampvis2, where data was first aggregated to family level and transformed to relative abundance.

Constrained ordination analysis was performed using two methods, canonical correspondence analysis and redundancy analysis (constrained principal component analysis), to investigate the relationships between sample types. Prior to the analysis ASVs < 0.1 relative abundance and samples with < 1000 sequences were removed. Data was first transformed using the Hellinger transformation and distance was measured using the Bray-Curtis method. Beta diversity statistical analysis was performed using analysis of variance (ANOVA) and permutational multivariate analysis of variance (PERMANOVA) (999 permutations) via the adonis function in the vegan package (Oksanen et al. 2020), (Hellinger transformation and Bray-Curtis distance measure), the F statistic and P value are presented (full statistical output available in Supplementary Material A.2.1). Hierarchical cluster analysis was performed using euclidean distance measure (average), where data was first transformed using Hellinger transformation.

The composition for taxa of interest was aggregated to family level and compared between the three sample types. Taxa of interest are defined as bacteria related to known pathogens associated with ticks and other vectors described globally, as outlined in Egan et al. (2020). Statistical analysis was performed using the Wilcoxon pairwise test (non-parametric). To identify shared ASVs (i.e. core microbiome) between sample types and species, data was transformed to relative abundance and assigned to best taxonomic hit using the microbiomeutilities R package. Similarities of the core microbiome between the three sample types were investigated using thresholds of abundance and prevalence between 0.0001–-0.001 and 0.001–-0.05 respectively. For subsequent comparisons of the core microbiome between host and tick species, abundance and prevalence levels were set at 0.001 and 0.05 respectively. Prevalence data for taxa of interest is presented for samples based on 16S rRNA metabarcoding. Where ticks have been pooled, prevalence is reported based on minimum infection rates (i.e. assumes only one positive tick in each pool) (Estrada-Pena et al. 2021).

Scripts for data analysis are available at siobhon-egan/wildlife-bacteria and siobhon-egan/wildlife-ticks and raw Illumina MiSeq data has been deposited in the European Nucleotide Archive under project accession number PRJEB46056 (16S rRNA bacteria) and PRJEB46056 (12S rRNA ticks).

3.3.5 Target PCRs

A nested PCR was used to amplify ~1.4 kb fragment of 16S rRNA gene for Neoehrlichia and Ehrlichia. Amplicon PCRs were carried out in 25 $\mu$L reactions each containing: 1X buffer (KAPA Biosystems, South Africa), 2.5 mM MgCl₂, 0.4 $\mu$M of each primer, 0.25 mM of each dNTP, 0.5 U of Taq (KAPA Biosystems, South Africa) and 2 $\mu$L of gDNA or 1 $\mu$L of primary product. Thermal cycling conditions were as follows; 95$^\circ$C for 3 mins, followed by 40 cycles of 95$^\circ$C for 30 secs, 48$^\circ$C (primary) or 54$^\circ$C (secondary) for 1 min, 72$^\circ$C for 2 mins; and a final extension of 72$^\circ$C for 5 mins.

For Bartonella, a PCR was used to amplify a ~370-450 bp fragment of the 16S rRNA - 23S internal transcribed spacer (ITS) region. Amplicon PCRs were carried out in 25 $\mu$L reactions each containing: 1X buffer (KAPA Biosystems, South Africa), 2.0 mM MgCl₂, 0.4 $\mu$M of each primer, 0.25 mM of each dNTP, 0.5 U of Taq (KAPA Biosystems, South Africa) and 2 $\mu$L of gDNA. Thermal cycling conditions were as follows; 95$^\circ$C for 3 mins, followed by 40 cycles of 95$^\circ$C for 15 secs, 66$^\circ$C 15 secs, 72$^\circ$C for 18 secs; and a final extension of 72$^\circ$C for 5 mins.

Amplicons were visualised on agarose gel and products of the expected size were excised with a sterile scalpel blade and purified using an in-house filtered pipette tip method (R. Yang et al. 2013). Sanger sequencing was performed at the Australian Genome Research Facility (Perth, Western Australia) on an Applied Biosystems 3730xl DNA Analyzer using BigDye(TM) Terminator v3.1 Cycle Sequencing Kit.

3.3.6 Phylogeny

Nucleotide sequences were inspected and quality filtered using Geneious 10.2.6 (https://www.geneious.com). Identity was confirmed using BLAST analysis (BLASTN 2.11.0+ (Zhang et al. 2000; Morgulis et al. 2008)) against NCBI nucleotide collection (nt) database. Sequences generated were aligned with references retrieved from GenBank (Benson et al. 2017) using MUSCLE (R. C. Edgar 2004) (Anaplasmataceae and Borrelia) or clustalW (Larkin et al. 2007) (Bartonella). The clustalW aligment method was identified as most suitable for Bartonella analysis due to variable length of the 16S rRNA and ITS regions targeted. Phylogenies were inferred using the maximum likelihood (ML) method. The optimal evolutionary model was selected using ModelFinder (Kalyaanamoorthy et al. 2017) based on the Bayesian information criterion. Phylogenetic analysis was performed in IQ-TREE v1.6.11 (Nguyen et al. 2015) and bootstrap support was calculated using ultrafast (UFBoot2) method with 10,000 replicates (Hoang et al. 2018).

3.4 Results

3.4.1 High-throughput sequencing

The hosts species sampled in the present study were the black rat (Rattus rattus), brown antechinus (Antechinus stuartii), brush-tailed possum (Trichosurus vulpecula), bush rat (Rattus fuscipes), chuditch (Dasyurus geoffroii), Rusa deer (Rusa timorensis), long-nosed bandicoot (Perameles nasuta), quenda (Isoodon fusciventer), rabbit (Oryctolagus cuniculus) and swamp rat (Rattus lutreolus) (Table 3.1). Ticks identified from wildlife hosts were Am. triguttatum, Ix. antechini, Ix. australiensis, Ix. holocyclus, Ix. tasmani and Ix. trichosuri. Molecular screening was performed on 257 ticks pooled into 205 gDNA extracts. Metabarcoding of ticks at the mitochondrial 12S rRNA gene confirmed morphological analysis and accurate differentiation of morphologically similar species. These results identified that mixed tick species were detected in five gDNA pools and in all cases, these were larvae of Ix. holocyclus and Ix. trichosuri (Table 3.3). A subset of representative 12S rRNA tick zOTU sequences were deposited under accession numbers MW665133–MW665150.

Table 3.3: Summary of tick species collected from wildlife. Sample numbers in table refer to number of genomic DNA extracts with total number of tick specimens included in parenthesis.
Host	Tick species	Larvae	Nymph	Male	Female
Black rat	Ix. holocyclus	6 (12)	5 (5)	0	0
	Ix. tasmani	30 (38)	31 (34)	0	0
	Ix. trichosuri	4 (9)	2 (3)	0	0
	Ix. holocyclus; Ix. trichosuri	1 (3)	0	0	0
Brown antechinus	Ix. antechini	0	0	0	1 (1)
	Ix. holocyclus	1 (1)	0	0	0
	Ix. tasmani	2 (2)	2 (2)	0	0
Brush-tailed possum	Am. triguttatum	2 (3)	0	0	0
	Ix. holocyclus	0	3 (3)	0	5 (5)
	Ix. tasmani	0	0	0	2 (2)
	Ix. trichosuri	1 (1)	2 (2)	0	18 (18)
Bush rat	Ix. tasmani	1 (1)	3 (3)	0	0
Chuditch	Am. triguttatum	1 (1)	0	0	0
Long-nosed bandicoot	Ix. holocyclus	8 (9)	5 (5)	5 (5)	20 (20)
	Ix. tasmani	5 (7)	12 (12)	0	0
	Ix. trichosuri	3 (3)	2 (2)	0	0
Quenda	Am. triguttatum	0	1 (1)	0	0
	Ix. australiensis	0	1 (1)	0	0
Rabbit	Ix. holocyclus	3 (6)	6 (10)	0	0
	Ix. tasmani	0	2 (2)	0	0
	Ix. trichosuri	1 (4)	4 (6)	0	0
	Ix. holocyclus; Ix. trichosuri	3 (12)	1 (3)	0	0

3.4.2 Bacteria 16S rRNA metabarcoding

A total of 536 gDNA extracts from 203 individual animals (159 blood, 205 tick pools, and 172 tissue samples) underwent bacterial profiling. An additional 65 control samples (EXB and NTC) were also sequenced to obtain background laboratory and reagent profiles. Bacterial 16S metabarcoding produced 102,503,643 raw sequences (Table A.4 and Figure A.2; samples = 96,969,355, controls = 5,534,288). Blood samples produced 35,160,160 sequences (mean = 185,053), tick samples produced 24,844,688 sequences (mean = 95,556) and tissue samples produced 36,964,507 sequences (mean = 197,671). Following quality filtering, denoising and merging a total of 64,186,307 sequences (66.2%) were obtained from the samples; mean for blood samples 121,376 (65.6%) mean for tick samples 70,650 (73.9%) and mean for tissue samples 121,689 (61.6%).

After contaminating taxa and controls were removed, a total of 31,194 bacterial ASVs were identified and, 5,243 ASVs had an abundance of over 1000 sequences. Rarefaction curves (Figure A.3) showed the number of ASVs generally plateaued at a depth of 60,000, 40,000 and 50,000 sequences for blood, tick and tissue samples respectively. Alpha diversity analysis (Figure 3.2) identified significant variation between sample types (Wilcoxon pairwise test, P < 0.001). The highest number of observed ASVs was in tissue samples (21,811), followed by blood (9,401) and ticks (3,568). Tissue samples showed the highest level of alpha diversity across all measures of alpha diversity. Ixodes australiensis and Ix. trichosuri had the highest alpha diversity of tick species sampled (Figure A.4). Once identified, contaminant taxa and controls were removed and the dominant phyla identified were Proteobacteria (37.3%), Firmicutes (27.8%), Actinobacteria (14.7%) and Bacteroidetes (10.7%) (Figures A.5 and A.6).

Boxplot of Alpha-diversity indices. Diversity indexes (A) Observed number of ASVs, (B) Chao1 index, (C) Shannon index and (D) inverse Simpson index. Boxplots and violin plots represent the distribution of diversity among samples within their category - blood (159), tick (205 pools) and tissue (172). Statistical analysis between sample types calculated using Wilcoxon pairwise (non-parametric) test with significance values indicated as follows: NS for P > 0.05; * for P <= 0.05; ** for P <= 0.01; *** for P <= 0.001; **** for P <= 0.0001.

Figure 3.2: Boxplot of Alpha-diversity indices. Diversity indexes (A) Observed number of ASVs, (B) Chao1 index, (C) Shannon index and (D) inverse Simpson index. Boxplots and violin plots represent the distribution of diversity among samples within their category - blood (159), tick (205 pools) and tissue (172). Statistical analysis between sample types calculated using Wilcoxon pairwise (non-parametric) test with significance values indicated as follows: NS for P > 0.05; * for P <= 0.05; ** for P <= 0.01; *** for P <= 0.001; **** for P <= 0.0001.

In blood samples, the most abundant families were Bacillaceae (22.9%), Bartonellaceae (11.3%), Anaplasmataceae (7.3%) and Mycobacteriaceae (2.6%) (Figure A.7). In tick samples the most abundant families were Coxiellaceae (33.8%), Midichloriaceae (31.0%), Mycobacteriaceae (6.0%), Flavobacteriaceae (2.6%) and Anaplasmataceae (1.8%) (Figure A.8). In tissue samples the most abundant families were Staphylococcaceae (8.2%), Mycobacteriaceae (6.5%), Ruminococcaceae (5.5%), Flavobacteriaceae (4.7%) and Anaplasmataceae (2.7%) (Figure A.9). Comparison of family taxa showed differences in abundance among sample types and species (and tick instar) (Figures 3.3 and 3.4).

Heatmap of the top 40 most prevalent bacterial family taxa identified in wildlife blood, tick and tissue samples. Data first transformed to relative sequence abundance.

Figure 3.3: Heatmap of the top 40 most prevalent bacterial family taxa identified in wildlife blood, tick and tissue samples. Data first transformed to relative sequence abundance.

Heatmap of the top 15 most prevalent bacterial family taxa identified in tick samples, showing abundance in species and instar. Data first transformed to relative sequence abundance. Tick samples consisting of mixed species pools were excluded.

Figure 3.4: Heatmap of the top 15 most prevalent bacterial family taxa identified in tick samples, showing abundance in species and instar. Data first transformed to relative sequence abundance. Tick samples consisting of mixed species pools were excluded.

The overlap of samples at the ASV level was investigated using various thresholds to define the ‘core microbiome’, with two out of three models showing only seven shared ASVs between all three sample types (Figure 3.5). Shared core microbiome ASVs of blood and tissue samples were compared between the four most abundant hosts samples and tick species (Figure 3.5).

$Venn diagram of core taxa present within samples presented as detection of amplicon sequence variants (ASVs) using abundance and prevalence. Data first transformed to relative abundance and ASVs assigned to best hit. (A) Impact of similarity based on detection thresholds of abundance and prevalence as (i) 0.001 and 0.05, (ii) 0.001 and 0.001, (iii) 0.0001 and 0.05. (B) Core taxa present within sample types using detection levels set at 0.001 and 0.05 for abundance and prevalence respectively for (i) blood (ii) tissue and (iii) ticks. For within sample comparison only the four most abundant vertebrate host/tick species were selected. Tick species abbreviations: \textit{Amblyomma triguttatum} (\textit{Am. tri}), \textit{Ixodes holocyclus} (\textit{Ix. hol}), \textit{Ixodes tasmani} (\textit{Ix. tas}), \textit{Ixodes trichosuri} (\textit{Ix. tri}). Tick samples consisting of mixed species pools were excluded.$

Figure 3.5: Venn diagram of core taxa present within samples presented as detection of amplicon sequence variants (ASVs) using abundance and prevalence. Data first transformed to relative abundance and ASVs assigned to best hit. (A) Impact of similarity based on detection thresholds of abundance and prevalence as (i) 0.001 and 0.05, (ii) 0.001 and 0.001, (iii) 0.0001 and 0.05. (B) Core taxa present within sample types using detection levels set at 0.001 and 0.05 for abundance and prevalence respectively for (i) blood (ii) tissue and (iii) ticks. For within sample comparison only the four most abundant vertebrate host/tick species were selected. Tick species abbreviations: (), (), (), (). Tick samples consisting of mixed species pools were excluded.

Ordination analysis revealed a distinct difference between the bacterial composition of blood, tick and tissue samples (Figure 3.6). Hierarchical cluster analysis (Figure A.10) and statistical analysis (PERMANOVA F = 36.209, P = 0.001, see Supplementary Material A.2.1) also supported this finding showing that there was a statistical difference in sample types, and there was little overlap in the microbiome composition of blood and tick samples. Tissue samples were identified as an intermediate sample type and showed some (albeit low) similarities to blood and tick samples. Ordination analysis of tick samples showed distinction between tick species (Figure 3.6). Tick species Ix. holocyclus and Ix. trichosuri showed the most similar bacterial composition with ordination analysis demonstrating overlap between these two species (PERMANOVA, F = 14.01, P = 0.001, see Supplementary Material A.2.1). Although sampling numbers were uneven among host species, which limited statistical inference, hierarchical cluster analysis did show that samples grouped were based on host (and tick) species (Figure A.11).

$Constrained ordination plot of microbiome composition from using constrained methods (A) Canonical Correspondence Analysis and (B) Redundancy Analysis (constrained Principal Component Analysis). (i) Ordination analysis of all samples N = 536 (159 blood, 205 tick pools, and 172 tissue) by sample type. (ii) Ordination of all tick samples by tick species (iii) Ordination of three co-habiting tick species in Sydney Northern Beaches area. Prior to the analysis ASVs < 0.1 relative abundance were removed. Analysis using the Hellinger transformation and Bray-Curtis distance measure. The relative contribution (eigen value) of each axis to the total inertia in the data as well as to the constrained space only, respectively, are indicated in percent at the axis titles. Tick species abbreviations: \textit{Amblyomma triguttatum} (\textit{Am. tri}),\textit{Ixodes antechinus} (\textit{Ix. ant}), \textit{Ixodes australiensis} (\textit{Ix. aus}), \textit{Ixodes holocyclus} (\textit{Ix. hol}), \textit{Ixodes tasmani} (\textit{Ix. tas}), \textit{Ixodes trichosuri} (\textit{Ix. tri}).$

Figure 3.6: Constrained ordination plot of microbiome composition from using constrained methods (A) Canonical Correspondence Analysis and (B) Redundancy Analysis (constrained Principal Component Analysis). (i) Ordination analysis of all samples N = 536 (159 blood, 205 tick pools, and 172 tissue) by sample type. (ii) Ordination of all tick samples by tick species (iii) Ordination of three co-habiting tick species in Sydney Northern Beaches area. Prior to the analysis ASVs < 0.1 relative abundance were removed. Analysis using the Hellinger transformation and Bray-Curtis distance measure. The relative contribution (eigen value) of each axis to the total inertia in the data as well as to the constrained space only, respectively, are indicated in percent at the axis titles. Tick species abbreviations: (), (), (), (), (), ().

3.4.3 Taxa of interest

Nine family taxa were chosen for comparison among sample types (Figure 3.7). Anaplasmataceae bacteria were most prevalent in blood samples (both by relative abundance and prevalence). Blood and tissue samples showed no statistical difference in relative abundance (Wilcoxon pairwise test, P = 0.26), with significantly lower abundance in tick samples (Wilcoxon pairwise test, P < 0.001). There was weak support for differences in Mycoplasmataceae between sample types (Wilcoxon pairwise test, P = 0.023–0.43). Midichloriaceae, Coxielleaceae, and Rickettsiaceae were all significantly more abundant in ticks than blood or tissue samples. Statistical comparisons of Borreliaceae and Francisellaceae were not made due to the small number of positive samples.

For subsequent prevalence analysis only taxa of interest that could be accurately defined to the genus level were retained (as per BLAST results). This identified 235 ASVs from 364 samples and 159 individuals see Table 3.4 and full sequence data and BLAST results in Supplementary Table E1.1.

Relative abundance for select bacteria taxa of interest (aggregated to family level) from wildlife samples. Statistical analysis between sample types calculated using Wilcoxon pairwise (non-parametric) test with significance values indicated as follows: NS for P > 0.05; * for P <= 0.05; ** for P <= 0.01; *** for P <= 0.001; **** for P <= 0.0001.

Figure 3.7: Relative abundance for select bacteria taxa of interest (aggregated to family level) from wildlife samples. Statistical analysis between sample types calculated using Wilcoxon pairwise (non-parametric) test with significance values indicated as follows: NS for P > 0.05; * for P <= 0.05; ** for P <= 0.01; *** for P <= 0.001; **** for P <= 0.0001.

Table 3.4: Presence of ‘taxa of interest’(a) by bacteria metabarcoding. Prevalence reported as proportion of individual hosts positive and proportion of samples positive in parenthesis (accounting for instances where more than one sample type was collected from an individual).
Taxa of interest	Sample type	Black rat	Brown antechinus	Brush tail possum	Bush rat	Chuditch	Deer	Long-nosed bandicoot	Quenda	Rabbit	Swamp rat
Anaplasmataceae	Blood	4/68 (4/84)	0/0 (0/0)	5/18 (6/21)	0/2 (0/2)	7/22 (9/24)	0/0 (0/0)	18/20 (22/25)	0/1 (0/2)	0/0 (0/0)	0/1 (0/1)
	Tick	3/54 (3/79)	1/4 (1/6)	4/15 (5/33)	0/3 (0/4)	0/1 (0/1)	0/0 (0/0)	8/31 (9/60)	0/2 (0/2)	0/7 (0/20)	0/0 (0/0)
	Tissue	6/71 (8/80)	0/2 (0/2)	4/16 (5/20)	0/3 (0/3)	1/9 (1/9)	3/3 (3/6)	18/40 (19/44)	0/0 (0/0)	0/7 (0/7)	0/1 (0/1)
Bartonella	Blood	22/68 (27/84)	0/0 (0/0)	7/18 (8/21)	2/2 (2/2)	0/22 (0/24)	0/0 (0/0)	4/20 (4/25)	0/1 (0/2)	0/0 (0/0)	0/1 (0/1)
	Tick	3/54 (4/79)	0/4 (0/6)	0/15 (0/33)	1/3 (1/4)	0/1 (0/1)	0/0 (0/0)	1/31 (1/60)	0/2 (0/2)	0/7 (0/20)	0/0 (0/0)
	Tissue	4/71 (5/80)	1/2 (1/2)	9/16 (11/20)	2/3 (2/3)	0/9 (0/9)	0/3 (0/6)	3/40 (3/44)	0/0 (0/0)	0/7 (0/7)	0/1 (0/1)
Borrelia	Blood	0/68 (0/84)	0/0 (0/0)	0/18 (0/21)	0/2 (0/2)	0/22 (0/24)	0/0 (0/0)	0/20 (0/25)	0/1 (0/2)	0/0 (0/0)	0/1 (0/1)
	Tick	0/54 (0/79)	0/4 (0/6)	0/15 (0/33)	0/3 (0/4)	0/1 (0/1)	0/0 (0/0)	0/31 (0/60)	0/2 (0/2)	0/7 (0/20)	0/0 (0/0)
	Tissue	7/71 (7/80)	0/2 (0/2)	0/16 (0/20)	1/3 (1/3)	0/9 (0/9)	0/3 (0/6)	0/40 (0/44)	0/0 (0/0)	0/7 (0/7)	1/1 (1/1)
Coxiellaceae	Blood	7/68 (7/84)	0/0 (0/0)	1/18 (1/21)	0/2 (0/2)	0/22 (0/24)	0/0 (0/0)	3/20 (5/25)	0/1 (0/2)	0/0 (0/0)	0/1 (0/1)
	Tick	50/54 (63/79)	3/4 (4/6)	8/15 (9/33)	3/3 (4/4)	0/1 (0/1)	0/0 (0/0)	22/31 (24/60)	0/2 (0/2)	2/7 (2/20)	0/0 (0/0)
	Tissue	33/71 (38/80)	0/2 (0/2)	2/16 (2/20)	2/3 (2/3)	0/9 (0/9)	0/3 (0/6)	20/40 (22/44)	0/0 (0/0)	0/7 (0/7)	0/1 (0/1)
Francisellaceae	Blood	0/68 (0/84)	0/0 (0/0)	0/18 (0/21)	0/2 (0/2)	0/22 (0/24)	0/0 (0/0)	0/20 (0/25)	0/1 (0/2)	0/0 (0/0)	0/1 (0/1)
	Tick	0/54 (0/79)	0/4 (0/6)	2/15 (2/33)	0/3 (0/4)	0/1 (0/1)	0/0 (0/0)	0/31 (0/60)	2/2 (2/2)	0/7 (0/20)	0/0 (0/0)
	Tissue	0/71 (0/80)	0/2 (0/2)	0/16 (0/20)	0/3 (0/3)	0/9 (0/9)	0/3 (0/6)	0/40 (0/44)	0/0 (0/0)	0/7 (0/7)	0/1 (0/1)
Midichloria	Blood	0/68 (0/84)	0/0 (0/0)	0/18 (0/21)	0/2 (0/2)	0/22 (0/24)	0/0 (0/0)	0/20 (0/25)	0/1 (0/2)	0/0 (0/0)	0/1 (0/1)
	Tick	17/54 (21/79)	2/4 (2/6)	12/15 (20/33)	0/3 (0/4)	0/1 (0/1)	0/0 (0/0)	19/31 (41/60)	0/2 (0/2)	5/7 (10/20)	0/0 (0/0)
	Tissue	5/71 (5/80)	0/2 (0/2)	1/16 (1/20)	0/3 (0/3)	0/9 (0/9)	0/3 (0/6)	3/40 (3/44)	0/0 (0/0)	0/7 (0/7)	0/1 (0/1)
Mycoplasma	Blood	1/68 (1/84)	0/0 (0/0)	0/18 (0/21)	0/2 (0/2)	0/22 (0/24)	0/0 (0/0)	0/20 (0/25)	0/1 (0/2)	0/0 (0/0)	1/1 (1/1)
	Tick	0/54 (0/79)	0/4 (0/6)	0/15 (0/33)	0/3 (0/4)	0/1 (0/1)	0/0 (0/0)	0/31 (0/60)	0/2 (0/2)	0/7 (0/20)	0/0 (0/0)
	Tissue	1/71 (1/80)	0/2 (0/2)	0/16 (0/20)	0/3 (0/3)	0/9 (0/9)	0/3 (0/6)	0/40 (0/44)	0/0 (0/0)	0/7 (0/7)	0/1 (0/1)
Rickettsia	Blood	0/68 (0/84)	0/0 (0/0)	0/18 (0/21)	0/2 (0/2)	0/22 (0/24)	0/0 (0/0)	0/20 (0/25)	0/1 (0/2)	0/0 (0/0)	0/1 (0/1)
	Tick	7/54 (7/79)	0/4 (0/6)	2/15 (2/33)	2/3 (3/4)	1/1 (1/1)	0/0 (0/0)	9/31 (10/60)	0/2 (0/2)	1/7 (1/20)	0/0 (0/0)
	Tissue	0/71 (0/80)	0/2 (0/2)	0/16 (0/20)	0/3 (0/3)	0/9 (0/9)	0/3 (0/6)	0/40 (0/44)	0/0 (0/0)	0/7 (0/7)	0/1 (0/1)
^a Taxa of interest are defined as bacteria related to known pathogens associated with ticks and other vectors described globally, as outlined in Egan et al. 2020.

3.4.3.1 Anaplasmataceae

The Anaplasmataceae family was identified in 95 samples from 59 individuals. The most prevalent genera were Neoehrlichia (43 individuals), followed by Ehrlichia (eight individuals) than Anaplasma (three individuals).

‘Candidatus Neoehrlichia arcana’ was identified in 52 samples: 29 blood, seven ticks (Ix. holocyclus), and 16 tissue. ‘Candidatus Neoehrlichia australis’ was identified in 33 samples: eight blood, six ticks (Ix. holocyclus, Ix. trichosuri), and 19 tissue. Sixteen samples had mixed infections of ‘Ca. N. arcana’ and ‘Ca. N. australis’ (six blood and ten tissue). No ticks had mixed Neoehrlichia infections. Three putative novel Anaplasmataceae species were identified from the present study, two Ehrlichia spp. and a single Neoehrlichia sp. The novel Neoehrlichia sp. was identified in three samples: a tissue sample from a black rat and a corresponding tick from the same individual (Ix. tasmani); and a tick from a brown antechinus (Ix. antechini). One novel Ehrlichia sp. was identified in five samples from two brush-tailed possums: two blood, one tick (Ix. holocyclus) and two tissue. A second novel Ehrlichia sp. was identified in seven samples from six chuditch (six blood and one tissue). Only one chuditch was positive in corresponding tissue and blood samples. Three deer tissue samples were infected with Anaplasma platys (Supplementary Table E1.1). Detection of Anaplasmataceae bacteria differed between tick species and instars. Ixodes holocyclus females had the highest prevalence of this taxon, followed by nymph and larval ticks, while no males were positive (Figure 3.4). Only Ix. tasmani nymph samples were positive and in the case of Ix. trichosuri, females and nymphs were identified with Anaplasmataceae and all larvae were negative.

Phylogenetic analysis was performed using a 1,244 bp fragment of the 16S rRNA locus (Figure 3.8). All ‘Ca. N. arcana’ amplified were identical to each other and showed 100% similarity to ‘Ca. Neoehrlichia arcana’ isolate HT94 from Ix. holocyclus, Australia (KU865447). Ehrlichia sequences from chuditch blood samples were 100% identical to each other and most similar to Ehrlichia sp. Anan from Ix. ovatus, Japan (AB028319, 98.94% similarity). The nearest named species was Ehrlichia muris (CP006917, 98.86% similarity) isolated from a wild mouse (Eothenomys kageus) in Japan (genetic distances Supplementary Table E1.2).

$Maximum likelihood (ML) phylogenetic reconstruction of \textit{Anaplasmataceae} based on a 1,244 bp alignment of the 16S rRNA locus. Substitution model K80 + F with 10,000 replicates. Node values correspond to bootstrap support where values > 0.7 are indicated by shaded circles. Number of substitutions per nucleotide position is represented by the scale bar. Sequences generated in the present study are in bold.$

Figure 3.8: Maximum likelihood (ML) phylogenetic reconstruction of based on a 1,244 bp alignment of the 16S rRNA locus. Substitution model K80 + F with 10,000 replicates. Node values correspond to bootstrap support where values > 0.7 are indicated by shaded circles. Number of substitutions per nucleotide position is represented by the scale bar. Sequences generated in the present study are in bold.

3.4.3.2 Midichloria

A total of 103 samples were positive for Midichloria which were identified in 94 tick samples from 62 individual animals. Positive tick species included Ix. holocyclus, Ix. tasmani, Ix. trichosuri from black rats, brown antechinuses, brush-tailed possums, long-nosed bandicoots and rabbits from New South Wales. While Midichloria was identified in 13 Ix. tasmani, the relative abundance was low (108–363 sequences; 0.1-0.4% relative abundance). Nine tissue samples (nine individuals) were positive for Midichloria however the number of sequences was generally low (111–14,002) when compared to tick samples; hosts included black rats (five), a brush-tailed possum (one) and long-nosed bandicoots (two). Only two individuals had positive tissue and tick samples, one brush-tailed possum and one long-nosed bandicoot. Midichloria was identified in all instars of Ix. holocyclus but only female and nymph stages of Ix. trichosuri (Figure 3.4). No blood samples were positive for Midichloria (Figure 3.3).

3.4.3.3 Coxiellaceae

The Coxiellaceae family was identified in 183 samples from 113 individuals. Eleven individuals had positive blood samples, although there was a relatively low number of sequences (100–2,819), from hosts including the black rat (seven), brush-tailed possum (one) and long-nosed bandicoot (three). One hundred and six ticks from 88 individuals were positive; Am triguttatum (one), Ix. tasmani (90) Ix. holocyclus (nine), and Ix. trichosuri. Fifty-seven individuals were positive for Coxiellaceae in tissue samples, including the black rat (33), brush-tailed possum (two), bush rat (two) and long-nosed bandicoot (20). Forty-five individuals had at least two sample types positive for Coxiellaceae, and only two black rats had positive blood, tick and tissue samples. All instars of Ix. tasmani sampled (female, nymph, larva) were positive for Coxiellaceae, but only female Ix. holocyclus were positive.

3.4.3.4 Rickettsiaceae

Rickettsia was identified exclusively in 24 tick samples from 22 individuals. Positive samples from Western Australia were all identified as originating in Am. triguttatum (three) from brush-tailed possums and quenda. Ticks identified as positive from New South Wales were: Ix. holocyclus (nine) from black rats and long-nosed bandicoots, Ix. tasmani (15) from black rats, bush rats, long-nosed bandicoots and rabbits, and Ix. trichosuri (three) from long nosed bandicoots. No tissue or blood samples were positive for Rickettsia.

3.4.3.5 Borrelia

A novel Borrelia sp. was identified in nine tissue samples: seven black rats, one bush rat and one swamp rat. The identity of these Borrelia sequences showed they were most similar to Borrelia R57 (AY626138, 95.68–96.70% similarity, Supplementary Table E1.1). Corresponding blood and ticks from these individuals were negative and no other Borrelia sequences were identified from any other samples.

Attempts to amplify and perform Sanger sequencing on the flaB gene using pan-Borrelia primers (S. M. Loh et al. 2016) were unsuccessful. Given these failed attempts and the small volume of gDNA from samples (tissue biopsy punch only 2 mm), an alternative approach was employed to amplify additional fragments of the 16S rRNA v3-4 hypervariable region on the Illumina MiSeq.

After bioinformatic analysis outlined above, a single ASV was generated from all nine samples. Phylogeny was reconstructed using the maximum likelihood method (GTR + R model) based on 426 bp alignment of the 16S rRNA locus (hypervariable region v3-4) (Figure 3.9). The novel Borrelia sequence showed 98.4% similarity with Borrelia sp. ALEPB216 (KF957671) from R. rattus in California, U.S.A (genetic distances Supplementary Table E1.3). The next most similar sequences were Borrelia sp. CA682 (KF957670, similarity 98.1%) and Borrelia sp. R57 (AY626138, similarity 97.6%). The topology showed these sequences formed a distinct ‘rodent’ clade, basal to the three major groups currently described (Figure 3.9). The nearest named species was Borrelia theileri (U38375, similarity 94.50%).

$Maximum likelihood (ML) phylogenetic reconstruction of \textit{Borrelia} genus based on a 431 bp alignment of the 16S rRNA locus (hypervariable region 3-4). Substitution model K2P + G4 with 10,000 replicates. Node values correspond to bootstrap support where values > 0.7 indicated by shaded circles. Number of substitutions per nucleotide position is represented by the scale bar. Sequence generated in the present study in bold.$

Figure 3.9: Maximum likelihood (ML) phylogenetic reconstruction of genus based on a 431 bp alignment of the 16S rRNA locus (hypervariable region 3-4). Substitution model K2P + G4 with 10,000 replicates. Node values correspond to bootstrap support where values > 0.7 indicated by shaded circles. Number of substitutions per nucleotide position is represented by the scale bar. Sequence generated in the present study in bold.

3.4.3.6 Bartonella

Bartonella sp. was identified in 69 samples from 47 individuals. Blood samples represented the majority of Bartonella positives with 35 individuals which included black rats (22), brush-tailed possums (seven), bush rats (two), and long-nosed bandicoots (four). Nineteen individuals had positive tissue samples; black rats (four), brown antechinus (one), brush-tailed possums (nine), bush rats (two) and long-nosed bandicoots (three); seven individuals had positive tissue and blood samples. Six tick samples were positive for Bartonella (from five individuals) with three tick species (Ix. holocyclus, Ix. tasmani and Ix. trichosuri) from black rats and long-nosed bandicoots. All tick samples that were positive also had corresponding positive blood or tissue samples.

Targeted sequencing on a subset of positive samples showed two genotypes of Bartonella queenslandensis that were 97.2% similar to each other with ten SNPs (genetic distances Supplementary Table E1.4. The next most similar sequence to genotype BR025 and BR048 was Ba. queenslandensis strain AUST/NH8 (EU111767) at 96.1% and 94.8% similarity respectively. Comparison of Ba. queenslandensis type sequences showed divergence up to 7.2%% (between EU111765 and EU111766) within the species. Five genotypes of Ba. coopersplainsensis were identified with 94.4–99.2% similarity to each other; and between 96.7–98.7% similarity to Ba. coopersplainsensis type sequence (EU111770) (Figure 3.10).

$Maximum likelihood (ML) phylogenetic reconstruction of \textit{Bartonella} genus based on a 533 bp alignment of the 16S-23S internal transcribed space (ITS) region. Substitution model HKY + F with 10,000 replicates. Node values correspond to bootstrap support where values > 0.7 indicated by shaded circles. Number of substitutions per nucleotide position is represented by the scale bar. Sequences generated in the present study in bold.$

Figure 3.10: Maximum likelihood (ML) phylogenetic reconstruction of genus based on a 533 bp alignment of the 16S-23S internal transcribed space (ITS) region. Substitution model HKY + F with 10,000 replicates. Node values correspond to bootstrap support where values > 0.7 indicated by shaded circles. Number of substitutions per nucleotide position is represented by the scale bar. Sequences generated in the present study in bold.

3.4.3.7 Francisella and Mycoplasma

Francisella was identified in four ticks, Am. triguttatum (three) and Ix. australiensis (one) from brush-tailed possums and quenda in Western Australia (Table 3.4, Supplementary Table E1.1). Mycoplasma was present in three samples, from three individuals: two blood samples from a black rat and a swamp rat and one tissue sample from a black rat (Table 3.4, Supplementary Table E1.1).

3.5 Discussion

This body of work represents the first large scale investigation into tick (and vector)-associated bacteria concurrently with their vertebrate hosts in urban and peri-urban locations, and identified novel tick-associated bacteria in Australian wildlife. Bacterial amplicon sequencing showed that there were limited taxa shared between blood, tick and tissue samples. The identification of bacterial taxa of interest was generally confined to certain hosts and/or sample types. To the best of the authors’ knowledge, this study represents the first report of Neoehrlichia and Ehrlichia from Australian wildlife (blood and/or tissue samples). The first identification of the rodent-associated Borrelia clade in Australia is also described. These findings show that geographic boundaries and host associations influence the prevalence of tick-associated bacteria. Importantly, these findings also provide support for the absence of northern hemisphere tick-borne pathogens, adding to the overwhelming body of evidence that any potential human tick-borne pathogens in Australia are likely endemic and genetically distinct from currently described pathogens in parts of North America and Europe.

Although a comparatively deep level of sequencing was achieved compared to similar studies (Rynkiewicz et al. 2015; Ge et al. 2018), rarefaction plots showed that blood samples did not reach the same asymptote as tick and tissue, indicating that deeper sequencing of these samples is needed to ensure complete coverage of bacterial diversity. Across all four alpha diversity measures, tissue samples showed the highest diversity and variance. In comparison, tick samples showed a less diverse bacterial composition. This is likely due to the dominance of endosymbiont bacteria in the tick microbiome, which has been demonstrated in previous studies of the same tick species as identified in the present study (Gofton, Oskam, et al. 2015; Siobhon L. Egan et al. 2020), and in tick species globally (Lalzar et al. 2012; J. C. Gil et al. 2020). The overall lack of similarity in bacterial communities between sample types observed in the present work has been reported in previous studies on small mammals and their vectors (ticks and fleas), indicating the complexity of microbial assemblages inhabiting these biological niches (Cohen et al. 2015; Rynkiewicz et al. 2015; H. Li, Li, and Qu 2018). The diagnostic value of 16S rRNA amplicon sequencing methods has not been thoroughly validated for tick-borne pathogens. One study however showed that ticks which tested positive for Borrelia using qPCR methods also had a high abundance of Borrelia using amplicon bacterial metabarcoding (J. L. H. Sperling et al. 2020). Therefore, even though Ix. holocyclus and Ix. tasmani had abundant endosymbionts, and the presence of bacterial species in lower abundance cannot be ruled out completely, it is likely that our finding of the absence of some bacterial species here reflects the true situation.

The microbiome of Ix. holocyclus and Ix. tasmani was dominated by a limited number of highly abundant ASVs. Other tick species, such as Ix. trichosuri and Ix. australiensis, showed a more diverse bacterial composition. The interaction between bacterial endosymbionts and pathogens in ticks is not yet fully understood. The presence of endosymbiont bacteria is thought to inhibit the colonisation by other potentially pathogenic bacterial species. For example, colonisation by the rickettsial endosymbiont (Rickettsia peacockii) in Dermacentor andersoni tick ovaries is thought to render the ticks more resistant to co-infection by a second rickettsial species (Rickettsia rickettsii), which is the main cause of Rocky Mountain Spotted Fever (W. Burgdorfer, Hayes, and Mavros 1981) (W. Burgdorfer, Hayes, and Mavros 1981; Macaluso et al. 2002).

Likewise, Dermacentor variabilis ticks infected with one strain of Rickettsia were shown to be less susceptible to secondary rickettsial infections (Macaluso et al. 2002). However, this theory is not uniformly accepted, and there is some rejection of Burgdorfer’s (W. Burgdorfer, Hayes, and Mavros 1981) initial findings. Instead, it is hypothesised that differences in R. peacockii prevalence in De. andersoni are due to microclimatic factors (Sam R. Telford 2009). Additionally, in recent years microbiome profiling has led to conflicting conclusions about the interaction between endosymbiont and pathogens (Sarah Irene Bonnet and Pollet 2020). Based on what is known from the northern hemisphere examples, tick species without a dominant endosymbiont may pose a greater risk of harbouring potentially pathogenic microbes (Sarah Irene Bonnet and Pollet 2020).

Equivalent studies from the northern hemisphere using 16S rRNA metabarcoding methods and wildlife surveillance to identify tick-borne pathogens are limited. Instead, pathogen surveillance studies are generally targeted to well-known pathogens. For example, studies of B. burgdorferi sensu lato have shown infection is not detrimental to the health of natural rodent reservoirs (Voordouw, Lachish, and Dolan 2015), which have been shown to harbour a diverse range of Borrelia species (Bunikis and Barbour 2005). In the present study, chuditch and long-nosed bandicoots had a high abundance of Anaplasmataceae in blood samples, while brush-tailed possums had the highest abundance of this taxon in tissue samples. In addition, Anaplasmataceae bacteria were mostly absent from larval ticks.

This suggests that the Anaplasmataceae bacteria identified are likely maintained via horizontal transmission (i.e. between host and tick) and transstadial routes, which has been observed with related bacteria in companion animals (Fourie et al. 2013; Almazan et al. 2020). In contrast, the identification of Midichloria across all life stages of Ix. holocyclus shows evidence of vertical (transovarial) transmission. Both horizontal and vertical transmission have been identified previously in other tick species with Midichloria (Mariconti et al. 2012; Serra et al. 2018). The absence of Midichloria in vertebrate hosts (blood and skin) in the present study suggests the horizontal infection route is not a major contributor to the presence of this endosymbiont in Ix. holocyclus and is consistent with previous studies (Mukhacheva and Kovalev 2017).

3.5.1 Taxa of interest

3.5.1.1 Anaplasmataceae

The present study has provided the first identification of Neoehrlichia spp. from Australian wildlife. Two Neoehrlichia species, ‘Ca. N. arcana’ and ‘Ca. N. australis’, were recently characterised from Ix. holocyclus along the east coast of Australia at a prevalence of 11.3% (44/391) (Gofton et al. 2016). This finding is consistent with other species of Neoehrlichia globally, which have been described circulating in a range of mammals and associated ticks (Kawahara et al. 2004; Muller et al. 2018). In contrast to Gofton et al. (2016), the results presented here show that ‘Ca. N. arcana’ is more prevalent than ‘Ca. N. australis’ in small mammals, which is consistent with a recent study of Australian wildlife ticks (Siobhon L. Egan et al. 2020). The main reservoir host for ‘Ca. N australis’ is likely to be medium to large macropods (e.g. eastern grey kangaroo (Macropus giganteus) or swamp wallaby (Wallabia bicolor)), which were not sampled in the present study.

The identification of three putative novel Anaplasmataceae species further demonstrates the high diversity of Anaplasmataceae from Australian ticks and wildlife. This finding adds to recent discoveries identifying new species from Am. triguttatum (Gofton et al. 2017), Bothriocroton concolor, (Siew-May Loh 2018), Ix. holocyclus (Gofton et al. 2016), and from the platypus (Ornithorhynchus anatinus) and its associated ticks Ix. ornithorhynchi (Gofton et al. 2018). The novel Anaplasmataceae species were less readily identified in tick vectors compared to ‘Ca. N. arcana’ and ‘Ca. N. australis’; however more targeted studies are needed to understand their prevalence and distribution.

The first identification of Anaplasma platys (formerly Ehrlichia platys) in Australia was made in 2001 (G. K. Brown, Roberts, and Aitken 2001), and has since been described from dogs throughout the country (G. Brown et al. 2006; Hii et al. 2012). To date it has not been identified in any native Australian wildlife, although Anaplasma platys has been reported from deer overseas (Yu et al. 2020). To the best of the authors’ knowledge, this is the first report of A. platys from wild deer in Australia. The vector of A. platys overseas, Rhipicephalus sanguineus (Snellgrove et al. 2020), is likely also the vector in Australia. The brown dog tick, syn. Rh. sanguineus tropical lineage was recently named Rhipicephalus linnaei (Jan Slapeta, Chandra, and Halliday 2021). Importantly, Ehrlichia canis, which was first identified in dogs and brown dog ticks in northern Australia during 2020, was not identified in ticks, tissue or blood from wildlife in this study.

3.5.1.2 Midichloria

While ‘Ca. M. mitochondrii’ (CMm) is known to be highly prevalent and abundant in ticks identified in the present study, the present study did not employ a blocking primer as used in previous studies (Gofton, Oskam, et al. 2015; Siobhon L. Egan et al. 2020). The justification for this was a deliberate choice to include CMm in the ‘taxa of interest’ investigation, i.e. provide insight into the ecology of this microbe. The absence of CMm in both blood and tissue was surprising, especially given the high proportion and abundance of the bacteria in Ix. holocyclus in the present study. Studies in France have identified CMm from roe deer blood samples (Capreolus capreolus) using molecular and serology techniques (Serra et al. 2018); while that study reported a higher sensitivity for serology, molecular methods were also successful at identifying infection. In vitro studies of Ixodes ricinus and rabbits suggest that the endosymbiont replicates within the vertebrate host (Cafiso et al. 2019). The bacterial microbiome of many tick species is generally dominated by the presence of one (or a few) highly abundant endosymbiont microbes (Clayton et al. 2015; Gofton, Oskam, et al. 2015; Varela-Stokes et al. 2017). This suggests that rare, less abundant bacteria are not easily identified using standard 16S rRNA bacterial metabarcoding methods. While the possibility cannot be ruled out, given the high abundance of CMm in tick samples, it is likely that the methods employed here would have identified it in wildlife hosts (either blood or tissue) if present. Therefore, given the absence of CMm in the hosts sampled in the present study, vertical transmission may be the primary source of Midichloria infection for Ix. holocyclus. Alternatively, the reservoir host for CMm may not have been sampled.

3.5.1.3 Coxiellaceae

All Coxiellaceae ASVs identified were most similar to members of the Rickettsiella genus and dominated the microbiome of Ix. tasmani. Coxiella was not detected in the present study. The absence of Coxiella from host blood and tissue samples likely reflects a low abundance of this bacteria in the small mammal hosts sampled. Serological methods are more commonly used to detect exposure to C. burnetii (A. Cooper et al. 2012) but lack the ability to detect acute infections. A recent study using molecular techniques detected C. burnetii in kangaroo meat samples (Shapiro et al. 2020). This suggests that animals sampled in the present study were negative for an active infection with Coxiella.

3.5.1.4 Rickettsiaceae

While several forms of rickettsiosis are amongst the few recognised human tick-borne diseases in Australia, remarkably little is known about the sylvatic life cycle of these intracellular bacteria. At least five agents of human rickettsial disease have been described in Australia (Tadepalli et al. 2021), together with several additional species and genotypes with unknown pathogenicity (Vilcins, Old, and Deane 2008; Vilcins, Fournier, et al. 2009; Tadepalli et al. 2021). Bandicoots and other small mammals have previously been described as reservoirs for Queensland tick typhus (R. australis), and Ix. holocyclus and Ix. tasmani are known vectors (Campbell and Domrow 1974; Sexton et al. 1991). Although the relationship between bandicoots and rickettsial infections is largely accepted there has been little validation of this association since the early publications.

In the present study, the genus Rickettsia was predominately identified from ticks, with a low prevalence and abundance in four tissue samples and was absent from all blood samples. Areas sampled in Sydney’s northern beaches are endemic for Queensland tick typhus and spotted fever. Therefore, the absence of spotted fever or typhus group Rickettsia from mammal hosts in the present study was unexpected. Recent studies utilising qPCR have reported Rickettsia in 15.4% (23/149) Ix. holocyclus from NSW (Stephen Graves et al. 2016) and 6.4% (13/203) ticks collected from wildlife in north-east Queensland (Hussain-Yusuf et al. 2020). Additional reports also recognise the presence of spotted fever group Rickettsia in Australian ticks (Cook, Scott, and Campbell 1967; Campbell and Domrow 1974; S. Graves, Unsworth, and Stenos 2006).

The lack of Rickettsia identified in ticks may be due to the highly abundant endosymbionts Midichloria or Coxiellaceae in Ix. holocyclus and Ix. tasmani respectively. While there are limited comparable studies that have used bacterial metabarcoding on wildlife blood samples, previous research has shown that molecular methods successfully identify Rickettsia in reservoir hosts (Ndeereh et al. 2017; Martello et al. 2019; Chaorattanakawee et al. 2021). Flinders Island Spotted Fever, caused by R. honei is associated with reptiles and the southern reptile tick Bt. hydrosauri (Stenos et al. 2003). A closely related strain, named R. honei marmionii, has been identified from various regions in Australia and is associated with Haemaphysalis novaeguineae (Lane et al. 2005; N. Unsworth et al. 2007). Interestingly no Haemaphysalis ticks were identified from hosts in the present study. Further investigation into the sylvatic cycle of Rickettsia in Australia would benefit from the inclusion of samples from reptile hosts and vertebrate hosts of Haemaphysalis ticks.

3.5.1.5 Borrelia

The novel Borrelia identified here came exclusively from black rats and was found at two sites in Sydney’s Northern Beaches area. All samples that tested positive were tissue, and no Borrelia sequences were identified from corresponding blood or tick samples.

The Borrelia species identified from black rats formed a clade with other described Borrelia species identified from rodents. Borrelia sp. R57 was identified from bank voles (Myodes glareolus; formerly Clethrionomys glareolus) and wood mice (Apodemus sylvaticus) in Spain at a prevalence rate of 8.5–12.0% using PCR methods (H. Gil et al. 2005). That study identified serological cross-reactivity of this species with B. burgdorferi sensu lato using western blot testing. A follow-up study (in Spain) further identified the genotype from 24.5% (62/253) of rodent tissue samples (Barandika et al. 2007). Identification of similar sequences was subsequently made from black rats (R. rattus) in California, U.S.A., with a prevalence of 43.5% (10/23) in tissue (skin) biopsy samples (Fedorova et al. 2014). The absence of this Borrelia species from our tick samples is interesting and mirrors the previous studies in Europe and the U.S.A (H. Gil et al. 2005; Fedorova et al. 2014). Despite multiple attempts during this study, amplification at the flaB locus was unsuccessful, consistent with previous research of members within the clade (H. Gil et al. 2005).

The extensive geographical range of this clade, yet tight host associations identified to date (i.e. confined to the Order Rodentia), raises questions regarding its origin. Particularly interesting is the paucity of reports globally, despite small mammals being the most widely studied reservoir host of tick-borne diseases (Zikeli and Zohdy 2020). Our phylogenetic analysis shows strong support for these Borrelia sequences forming a distinct clade and is consistent with other studies (Margos et al. 2019). Its basal phylogenetic position may provide important insights into the evolutionary history of the Borrelia genus. The lack of its identification in tick samples means that the vector for this Borrelia sp. remains unknown, and the transmissibility and pathogenicity of this microbe requires further research. Importantly, no identifications were made of the recently characterised Borrelia species from Australia, ‘Ca. B. tachyglossi’ (S. M. Loh et al. 2017) or reptile Borrelia sp. (Panetta et al. 2017), nor any sequences from the B. burgdorferi sensu lato group.

Currently there is no evidence that the novel rodent-associated Borrelia genotype identified in the present study is transmissible to humans or capable of causing disease. Phylogenetic reconstruction indicates that it is distinct from the B. burgdorferi sensu lato and relapsing fever groups (Margos et al. 2019). The absence of detection from tick samples and its identification only in rodents, more likely implies another vector that presumably has a tight host relationship with Rattus species.

3.5.1.6 Bartonella

Bartonella spp. was identified in a relatively high proportion of wildlife blood samples, and most prevalent in blood samples from black rats (31.0%). Bartonella coopersplainsensis and Ba. queenslandensis were described from native Australian rodents (Melomys spp. and Rattus spp.) in Queensland (Gundi et al. 2009). Since then, a number of novel Bartonella strains have been identified from Australia, including novel species described from marsupials and their ectoparasites (ticks and fleas) (Vilcins, Kosoy, et al. 2009; Kaewmongkol, Kaewmongkol, Owen, et al. 2011; Kaewmongkol, Kaewmongkol, McInnes, et al. 2011). While generally considered to be transmitted by fleas, reports of Bartonella species isolated from ticks are increasing, as studies take a broader approach to surveillance of vector-associated microbes (Rynkiewicz et al. 2015; Asyikha, Sulaiman, and Mohd-Taib 2020; Krol et al. 2021). Identification of Bartonella in fed ticks does not mean they are a viable vector, but it does highlight that future studies should include ticks when exploring candidate vectors. The high identification of Bartonella in rodents has been reported in other metabarcoding studies from blood samples (Cohen et al. 2015).

3.5.1.7 Francisella and Mycoplasma

Francisella was identified exclusively from ticks in Western Australia, Am. triguttatum and Ix. australiensis. Previous bacterial profiling of Am. triguttatum has shown a high proportion of Francisella species (Gofton, Doggett, et al. 2015; Siobhon L. Egan et al. 2020). Given its high abundance in tick samples, the absence of Francisella in corresponding host may indicate a lack of transmission. Tularemia (caused by Francisella tularensis) can cause disease in humans and animals. While the infection is more commonly acquired via handling infected animal material or aerosol inhalation, F. tularensis has been reported in ticks (Regier et al. 2019). Francisella has been reported in several tick species (Garcia-Vozmediano et al. 2020). While the species identified in the present study is not closely related to F. tularensis, considering the significance of this genera (Regier et al. 2019), further characterisation should be of high importance.

The present study did not identify any haemoprotic Mycoplasma species. The Mycoplasma sequences identified were most closely related to Mycoplasma arthritidis, which is recognised as a common pathogen in rats (Dybvig et al. 2008). It was identified in two black rats and a single identification in the native swamp rat.

3.5.2 Limitations and recommendations

While the cost of high-throughput sequencing platforms continues to decline, analysis of large numbers of samples is still constrained by resources and funds. Amplicon high throughput sequencing remains a cost-effective approach for characterising bacterial microbiomes. With the advent of long-read sequencing technologies (e.g. PacBio), more studies are moving towards obtaining full-length 16S rRNA sequences which allow for more robust phylogenetic reconstruction and in many cases identification of bacteria to strain/genotype level (Earl et al. 2018; Callahan et al. 2019). In addition to describing the bacteria present, the idea of the functional tick microbiome has been suggested (Estrada-Pena, Cabezas-Cruz, and Obregon 2020). While adaption of these methods can seem promising, their use outside of human microbiome studies is questionable due to the reliability of databases (Sarah Irene Bonnet and Pollet 2020; Sun, Jones, and Fodor 2020). An alternative approach such as the “common core” or ecological “core microbiome” may be more suitable in the context of tick microbiome research (Risely 2020). Studies incorporating transcriptomics will likely also prove useful for understanding the functional components of the microbiome.

While the present study reported on the relative abundance of microbes between ticks sampled, it is well known that differences in the 16S copy number limit inferences (Brooks et al. 2015; Louca, Doebeli, and Parfrey 2018). In the context of the present study, relative abundance of taxa was compared among samples, as it was considered a suitable measure given samples were processed in the same manner. Calculations of absolute abundance using 16S rRNA metabarcoding are challenging. While copy number is one factor that can be modelled it relies on information that in many cases is not accurately known. Using mock communities can help assess bias without baseline data on copy number, however, these analysis methods are still limited.

The sensitivity of molecular diagnostic testing for tick-borne pathogens has been shown to vary greatly. Factors that can influence molecular detection of pathogens include level of parasitaemia/bacteremia, sample type, time of sample collection from initial infection, and and pre-analytical preparation techniques; additionally, this is compounded by differences between taxa due to host-pathogen interactions (Fuente et al. 2017; Estrada-Pena, Fuente, and Cabezas-Cruz 2017; Sarah I. Bonnet, Nijhof, and Fuente 2018). A strong advantage of using generic metabarcoding assays is their unbiased approach; in the context of the present study, the most important factor was to achieve unbiased characterisation of bacterial communities from Australian wildlife. However, it is acknowledged that the use of bacterial metabarcoding for the identification of low abundant and rare organisms is limited. The lack of recognised endemic bacterial pathogens such as Coxiella burnetii, and Rickettsia spp. was surprising. Future studies should use serological and culture methods in addition to these highly specific and sensitive molecular tools (e.g. qPCR).

Another limitation of the present study was the bias in animals sampled. Small mammals were chosen as it was hypothesised they would have the highest prevalence of potential tick-borne pathogens based on what is known from the northern hemisphere (Bastian et al. 2012; Bouchard et al. 2013; Swei et al. 2019). Further sampling with expanded host and spatial aspects are vital to understand the full epidemiology of tick (and vector) associated microbes. Likewise, alternative trapping methods (e.g., differing trap type, bait) may expand the suite of hosts captured, importantly those with territorial overlap (e.g., semi-arboreal marsupials).

In the present study, most taxa of interest were confined to specific sample types. For example, Borrelia was identified from tissue samples of rats but absent from any blood or tick samples. This finding has also been reported in white-footed mice (Peromyscus leucopus) (Zawada et al. 2020). While in many cases disease surveillance is targeted to known pathogens where sampling protocols and diagnostics are well defined, in scoping studies such as the present one, a range of sample types are essential. This is evident in the different ways pathogens act within their host. As the first line of defense, the innate immune system has a key role in protecting against cutaneous and systemic infection (Coates, Blanchard, and MacLeod 2018). Alternatively pathogens can quickly evade host innate immune responses by rapidly disseminating into the blood stream after invasion (D. Kumar et al. 2015; Hyde 2017). In the context of Australia however, many wildlife species are endangered and therefore their capture and sampling in the wild is subject to strict ethical conditions and government regulations. As such, collection of multiple sample types is restricted to those requiring only a minimally invasive method and is also dependent on the availability of trained personnel and resources. In many instances it is standard procedure to collect a tissue biopsy for animal genetic (population) purposes, this presents a simple alternative to other procedures which may require field anaesthesia. The results here highlight the value of this sample type for use in pathogen surveillance.

3.6 Conclusions

Increased land-use and climate change will continue to affect the spread of vectors and associated diseases globally. There is an urgent need for expanded disease surveillance, particularly at the urban interface with wildlife. Through collaborative approaches, discoveries of novel microbes and important epidemiological links will be critical to ensure adequate and timely responses to human and animal diseases. Investigations into taxa of interest in the present study showed Anaplasmataceae bacteria had a relatively high prevalence in the vertebrate hosts examined here and a wide geographical distribution. This group of bacteria is therefore a viable candidate for continued research into Australian zoonotic tick-borne diseases.

2 An update on records of Australian ticks

4 Haemoprotozoa surveillance