Vector-borne pathogens exist in obligate transmission cycles between vector and reservoir host species. Host and vector shifts can lead to geographic expansion of infectious agents and the emergence of new diseases in susceptible individuals. Three bacterial genospecies (Borrelia afzelii, Borrelia bavariensis, and Borrelia garinii) predominantly utilize two distinct tick species as vectors in Asia (Ixodes persulcatus) and Europe (Ixodes ricinus). Through these vectors, the bacteria can infect various vertebrate groups (e.g. rodents, birds) including humans where they cause Lyme borreliosis, the most common vector-borne disease in the Northern hemisphere. Yet, how and in which order the three Borrelia genospecies colonized each continent remains unclear including the evolutionary consequences of this geographic expansion. Here, by reconstructing the evolutionary history of 142 Eurasian isolates, we find evidence that the ancestors of each of the three genospecies likely have an Asian origin. Even so, each genospecies studied displayed a unique sub-structuring and evolutionary response to the colonization of Europe. The pattern of allele sharing between continents is consistent with the dispersal rate of the respective vertebrate hosts, supporting the concept that adaptation of Borrelia genospecies to the host is important for pathogen dispersal. Our results highlight that Eurasian Lyme borreliosis agents are all capable of geographic expansion with host association influencing their dispersal; further displaying the importance of host and vector association to the geographic expansion of vector-borne pathogens and potentially conditioning their capacity as emergent pathogens.

Isolates, culturing, and DNA extraction

The major aim of this study was to estimate the direction of colonization in each Borrelia genospecies between Asia and Europe. To achieve this though, we needed to create an isolate library of the three genospecies including isolates arising from wild-caught ticks, where possible. Therefore, this study utilized DNA of 136 Eurasian Borrelia isolates coming from B. afzelii (Asian, n=20; European, n=13), B. bavariensis (Asian, n =27; European, n=19), and B. garinii (Asian, n=25; European, n=32). Of these, 52 are novel Borrelia isolated from wild-caught ticks collected either in Japan (n=43) or Germany (n=9) in the years 2018–2019 (see Text S1 & Tables S2, S3, S4). Additionally, 55 European isolates (B. afzelii, n=11; B. garinii, n=25; and B. bavariensis, n=19) and 12 Japanese isolates (B. bavariensis, n=8; B. garinii, n=4) were provided by the German National Reference Center for Borrelia at the Bavarian Food and Health Safety Authority and the National Institute of Infectious Disease respectively. These additional isolates predominantly come from LB patients (n=59) with a few isolated from wild-caught ticks (n=8). Finally, Russian B. bavariensis (n=6) (Becker et al., 2020) and B. garinii (n=11) isolates arising from wild-caught ticks were included in the study (see Text S1). For three isolates (UO2, UO3, UO4), the source material is not known. For all information on isolates, including origin and source material, refer to Table S1.

Borrelia isolates were cultured either in inhouse-made MKP (Preac-Mursic et al., 1986) (all European isolates) or inhouse-made BSK-H (Pollack et al., 1993; Takano et al., 2014) (all Russian and Japanese isolates) medium according to standard procedures (Pollack et al., 1993; Preac-Mursic et al., 1986) until the cultures reached a density of at least 10⁸ cells per mL, at which point whole genomic DNA was extracted. Genomic DNA from all European isolates was extracted using a Maxwell^® 16 LED DNA kit (Promega, Madison, WI, USA) and from all Japanese and Russian isolates using the Wizard^® Genomic DNA Purification Kit (Promega, WI, USA). DNA quality (260/280) and concentration were measured using a NanoDrop^® 1000 photometer (Thermo Fisher Scientific, Waltham, MA, USA) and a Qubit^® 3.0 fluorometer (Thermo Fisher Scientific, MA, USA), respectively.

Previously isolated samples provided by collaboration partners were assigned to genospecies through amplifying and sequencing the eight MLST genes (Margos et al., 2008). Isolates produced during this study were assigned to genospecies through amplifying one MLST gene, recG, and by comparing this sequence to the PubMLST database for Borrelia (https://pubmlst.org/borrelia). The isolate was assigned to the closet allele type match. This approach has been used in previous work reliably to assign Borrelia genospecies (e.g. Rollins et al., 2021). For further details regarding isolation and processing see Text S1.

Whole genome sequencing, assembly, and determining orthologus genes

Within the course of this study, we aimed to reconstruct the evolutionary history of B. afzelii, B. bavariensis, and B. garinii for which we needed to focus on the portions of the genome which are stable and homologous across all isolates. For this, we aimed to reconstruct full chromosome sequences through whole genome sequencing for each isolate, as the chromosome represents a large proportion of the Borrelia genome and the majority of the sequence (>93%) is stable (Schwartz et al., 2021; Walter et al., 2017). For all samples, libraries were produced according to the Nextera XT sample preparation guide (Illumina, San Diego, CA, USA). Library quality was checked using an Agilent TapeStation 2200 (Agilent, Santa Clara, CA, USA) before being sequenced using an Illumina MiSeq platform according to standard protocol (Illumina, USA) that produced paired-end reads of 250bp. Illumina reads were first trimmed for Illumina MiSeq adapter sequences using Trimmomatic v. 0.38 (Bolger et al., 2014a, 2014b) before being assembled using SPAdes v. 3.13.0 (Bankevich et al., 2012), which has been shown to be the best option for de novo assemblies of Borrelia genomes (Becker et al., 2020). Pacific Bioscience sequences were obtained for three B. bavariensis isolates (PBi, A104S, and NT24) (Becker et al., 2020) and three B. garinii isolates (PHeI, PBr, and NT31; see Suppl. Met.). Additionally, three B. afzelii chromosomes were downloaded from GenBank for use as references and inclusion in all analyses: PKo (CP009058.1), K78 (CP002933.1), and ACA-1 (NZ_ABCU00000000.2). SPAdes contigs were then mapped to reference chromosomes using NUCmer v. 3.23 from the package MUMmer (Delcher et al., 2002; Kurtz et al., 2004). Final chromosomes were produced according to the mapping protocol outlined in Becker et al. (2020) (see Suppl. Met.). Three additional B. bavariensis chromosomes were downloaded from GenBank and used in further analyses: SZ (CP007564.1), BgVir (CP003151.1), and NWJW1 (CP003866.1).

Final chromosomes were uploaded and annotated using the RAST annotation server (Aziz et al., 2008; Overbeek et al., 2014) for which proposed coding sequences were extracted. Orthologous sequences were determined by using the CRBHits package (Ullrich, 2020) as implemented in R (R Core Team, 2019). Briefly, all coding sequences for each chromosome were compared pairwise to all other chromosomes using the crb2rbh function. This information was fed into the integrated DAGchainer (Haas et al., 2004) command to create links between syntenous genes. Synteny member groups were then determined using the “cluster_infomap” command as implemented from the igraph R package (Csárdi and Nepusz, 2006) on a matrix of gene names and DAGchainer links. FASTA files for all synteny member groups were generated and aligned using MUSCLE (Edgar, 2004a, 2004b) for all gene copies and each genospecies. For each alignment, nucleotide diversity (π) (Nei, 1987) was estimated in R (R Core Team, 2019) using the package pegas (Paradis, 2010) as well as the proportion of isolates carrying the synteny member group. 

Recombination analysis

Recombination is known to be low on the Borrelia chromosome (Gatzmann et al., 2015; Schwartz et al., 2021; Walter et al., 2017). Even so, recombinant areas of the genome are inherited through horizontal gene transfer (Arnold et al., 2022) and could therefore bias phylogenetic reconstruction. For this reason, we aimed to identify recombinant regions along the chromosome to remove prior to phylogenetic reconstruction but also to determine what genes may be influenced by recombination in our dataset. Final assembled chromosomes and GenBank sequences (n=142) were aligned using MAFFT v. 7.407 (Katoh et al., 2002; Katoh and Standley, 2013). To determine recombinant regions along the chromosomes, we applied the four-gamete condition (Hudson & Kaplan, 1985) to the full chromosome alignment, as described in Gatzmann et al. (2015).

The ordered list of segregating sites along the chromosome was divided into blocks containing the same number of SNPs (n = 12). Each pair of SNPs in each block was then assessed if the four-gamete condition was violated or not. The within-block average and standard deviation were then calculated and averaged across all blocks and used as a measure of background violation due to double hits or back mutations. To single out SNP blocks that were most likely under recombination, we then calculated all pairwise comparisons between blocks and recorded the violation score. This score was then averaged over all comparisons for a specific block. Blocks were considered recombinant if:

and non-recombinant if:

where x_i is the average violation of block i over all comparisons, μ_within is the average within-block violation score, and sd_within is the standard deviation of within-block violation score.

We ran the analysis seven times using each time an alignment containing one, two or all genospecies. In this way, all crosswise comparisons were accounted for. The proportion of recombinant blocks was calculated for each comparison as the number of blocks meeting the above requirement divided by the total number of SNP blocks for the alignment in question. Recombinant blocks were extracted and compared using BLAST v.2.8.1 (Altschul et al., 1990; Camacho et al., 2009) (algorithm: blastn) to the full chromosome for each isolate. BLAST start and end positions were then compared to the RAST annotation files to determine if the recombining region fell within a coding (synteny member group) or non-coding region of the chromosome. Synteny member groups identified as being influenced by recombinant SNP blocks in all crosswise comparisons were considered to be influenced by horizontal gene transfer and were analyzed phylogenetically in BEAST v2.6.6 (Bouckaert et al., 2019) with the following parameters: coalescent model with constant size, strict clock (Drummond et al., 2006) with a clock-rate fixed to 1, GTR substitution model with four gamma categories (Tavaré, 1986). Potential recombination events were checked by plotting individual gene phylogenies against the full chromosome tree produced without potential recombining regions (see following methods section), using the “cophylo” in the R package phytools (Revell, 2012).

Phylogenetic reconstruction including biogeographical inference

In the aim to estimate the pattern of Borrelia colonization across Eurasia, we needed to estimate the probability of each ancestor having a European or Asian origin. This required phylogenetic reconstruction considering biogeographical inference for which we could only use homologous, non-recombining portions of the chromosome. For this, chromosome regions identified as recombinant in the previously described analysis were removed and the sequences were realigned using MAFFT v. 7.407 (Katoh et al., 2002; Katoh and Standley, 2013) (final alignment length: 936,908 bp). Phylogeny reconstruction was performed in BEAST v2.6.6 (Bouckaert et al., 2019) with the following parameters: coalescent model with constant size, strict clock (Drummond et al., 2006) with a clock-rate fixed to 1, GTR substitution model with four gamma categories (Tavaré, 1986). Geographic location was modeled as a discrete trait (Europe or Asia) and included in the phylogenetic inference utilizing the default settings for a discrete character (i.e. symmetric mutation death model with one gamma category) (Lemey et al., 2009; Wallace et al., 2007). Three independent runs were performed each with 50 million steps chain with a relative burn-in of 20% before selecting the best tree with TreeAnnotator v. 1.10.4 (Drummond and Rambaut, 2007). Convergence of parameters was checked with Tracer v. 1.7.1 (Rambaut et al., 2018). Convergence to a single topology in all three independent runs was checked manually in FigTree v. 1.4.4 (http://tree.bio.ed.ac.uk/software/figtree/), which was also used to plot the tree shown in Figure 2.

Population genetic analyses

An additional aspect in terms of estimating the pattern of colonization between Asia and Europe was to test how genetic diversity within each genospecies varied between the continents. This was in the aim to better understand how proposed host and vector associations have potentially influence the evolution and dispersal of these bacteria. We further extend this to ask if the inclusion of non-random samples could bias population genetic statistics as our data set includes non-randomly sampled isolates (predominantly isolated from LB patients) but also randomly sampled isolates arising from wild-caught ticks. All statistical analysis was performed in R v. 3.6.1 (R Core Team, 2019). Genetic diversity (π) (Nei, 1987) and Tajima’s D test statistic (Tajima, 1989) were estimated in the package pegas (Paradis, 2010). Analysis of molecular variance (AMOVA) (Excoffier et al., 1992) was performed using the package poppr (Kamvar et al., 2015) whereas F_ST (Nei, 1987) and D_XY (Hudson, Slatkin, & Maddison, 1992) were estimated with the package PopGenome (Pfeifer et al., 2014).

Identification of plasmid content through plasmid partitioning genes

The primary aim of this study was to utilize full chromosome sequences to estimate the direction of Borrelia colonization between Europe and Asia as well as quantify the variability in genetic diversity between populations. Even so, the Borrelia genome contains many accessory plasmid types which have also been shown to be influenced by demographic events such as the bottleneck observed in B. bavariensis where the overall diversity of plasmid types also decreased due to colonizing a novel tick vector (Becker et al., 2020). To determine if this was also the case for the other two Eurasian genospecies, we estimated plasmid content based on the number of unique plasmid partitioning genes present in each assembly. These partitioning genes have been shown to be plasmid-specific and exist as single copies in Borrelia isolates (Casjens et al., 2012; Casjens and Huang, 1993; Fraser et al., 1997). Identification of plasmid partitioning genes was performed as outlined in Becker et al. (2020) (see Suppl. Met.). Briefly, we used BLAST v.2.8.1 (Altschul et al., 1990; Camacho et al., 2009) (algorithm: blastn) to search for the presence of plasmid partitioning genes of the PFam32, 49, 50, and 57.62 families in the assembled SPAdes contigs. Hits were removed if they did not cover more than half the length of the references and had lower than 80% percent identity. After curation, we defined a plasmid being present if at least one of the partitioning genes was present in the assembled contigs.

Standard two-side, unpaired t-tests were run on plasmid number between genospecies comparing the two geographic populations using the function t.test from the base R package (R Core Team, 2019). Classical multidimensional scaling (MDS) was run using the cmdscale function using the base R package on a distance matrix calculated from the binary presence/absence plasmid data per isolate. Further effects on plasmid content were tested using a generalized linear mixed effects model assuming a Poisson error distribution using the glmer function from the package lme4 (Bates et al., 2015). Fixed effects were included for sample origin (Asia vs. Europe) and source (human vs. tick isolate) and genospecies was fitted as a random effect. Mean estimates and their 95% credible intervals were estimated based on 5,000 simulations using the sim function from the package arm (Gelman and Su, 2016). Residual error was calculated according to Nakagawa & Schielzeth (2010).

References

Altschul, S.F., Gish, W., Miller, W., Myers, E.W., Lipman, D.J., 1990. Basic local alignment search tool. J. Mol. Biol. 215, 403–410. https://doi.org/10.1016/S0022-2836(05)80360-2

Arnold, B.J., Huang, I.T., Hanage, W.P., 2022. Horizontal gene transfer and adaptive evolution in bacteria. Nat. Rev. Microbiol. 20, 206–218. https://doi.org/10.1038/s41579-021-00650-4

Aziz, R.K., Bartels, D., Best, A.A., DeJongh, M., Disz, T., Edwards, R.A., Formsma, K., Gerdes, S., Glass, E.M., Kubal, M., Meyer, F., Olsen, G.J., Olson, R., Osterman, A.L., Overbeek, R.A., McNeil, L.K., Paarmann, D., Paczian, T., Parrello, B., Pusch, G.D., Reich, C., Stevens, R., Vassieva, O., Vonstein, V., Wilke, A., Zagnitko, O., 2008. The RAST Server: rapid annotations using subsystems technology. BMC Genomics 9, 75. https://doi.org/10.1186/1471-2164-9-75

Bates, D., Maechler, M., Bolker, B., Walker, S., 2015. Fitting Linear Mixed-Effects Models Using lme4. J. Stat. Softw. 67.

Becker, N.S., Rollins, R.E., Nosenko, K., Paulus, A., Martin, S., Krebs, S., Takano, A., Sato, K., Kovalev, S.Y., Kawabata, H., Fingerle, V., Margos, G., 2020. High conservation combined with high plasticity: Genomics and evolution of Borrelia bavariensis. BMC Genomics 21, 702. https://doi.org/10.1186/s12864-020-07054-3

Bolger, A.M., Lohse, M., Usadel, B., 2014a. Trimmomatic: A flexible trimmer for Illumina sequence data. Bioinformatics 30, 2114–2120. https://doi.org/10.1093/bioinformatics/btu170

Bolger, A.M., Lohse, M., Usadel, B., 2014b. Genome analysis Trimmomatic: a flexible trimmer for Illumina sequence data. Bioinformatics 30, 2114–2120. https://doi.org/10.1093/bioinformatics/btu170

Bouckaert, R., Vaughan, T.G., Barido-Sottani, J., Duchêne, S., Fourment, M., Gavryushkina, A., Heled, J., Jones, G., Kühnert, D., De Maio, N., Matschiner, M., Mendes, F.K., Müller, N.F., Ogilvie, H.A., du Plessis, L., Popinga, A., Rambaut, A., Rasmussen, D., Siveroni, I., Suchard, M.A., Wu, C.-H., Xie, D., Zhang, C., Stadler, T., Drummond, A.J., 2019. BEAST 2.5: An advanced software platform for Bayesian evolutionary analysis. PLOS Comput. Biol. 15, e1006650.

Camacho, C., Coulouris, G., Avagyan, V., Ma, N., Papadopoulos, J., Bealer, K., Madden, T.L., 2009. BLAST+: Architecture and applications. BMC Bioinformatics 10, 1–9. https://doi.org/10.1186/1471-2105-10-421

Casjens, S., Huang, W.M., 1993. Linear chromosomal physical and genetic map of Borrelia burgdorferi, the Lyme disease agent. Mol. Microbiol. 8, 967–980. https://doi.org/10.1111/j.1365-2958.1993.tb01641.x

Casjens, S., Mongodin, E.F., Qiu, W.G., Luft, B.J., Schutzer, S.E., Gilcrease, E.B., Huang, W.M., Vujadinovic, M., Aron, J.K., Vargas, L.C., Freeman, S., Radune, D., Weidman, J.F., Dimitrov, G.I., Khouri, H.M., Sosa, J.E., Halpin, R.A., Dunn, J.J., Fraser, C.M., 2012. Genome stability of Lyme disease spirochetes: Comparative genomics of Borrelia burgdorferi plasmids. PLoS One 7. https://doi.org/10.1371/journal.pone.0033280

Csárdi, G., Nepusz, T., 2006. The igraph software package for complex network research Title. InterJournal, Complex Syst. 1695.

Delcher, A.L., Phillippy, A., Carlton, J., Salzberg, S.L., 2002. Fast algorithms for large-scale genome alignment and comparison. Nucleic Acids Res. 30, 2478–2483. https://doi.org/10.1093/nar/30.11.2478

Drummond, A.J., Ho, S.Y.W., Phillips, M.J., Rambaut, A., 2006. Relaxed phylogenetics and dating with confidence. PLoS Biol. 4, 699–710. https://doi.org/10.1371/journal.pbio.0040088

Drummond, A.J., Rambaut, A., 2007. BEAST: Bayesian evolutionary analysis by sampling trees. BMC Evol. Biol. 7, 1–8. https://doi.org/10.1186/1471-2148-7-214

Edgar, R.C., 2004a. MUSCLE: a multiple sequence alignment method with reduced time and space complexity. BMC Bioinformatics 5, 113. https://doi.org/10.1186/1471-2105-5-113

Edgar, R.C., 2004b. MUSCLE: multiple sequence alignment with high accuracy and high throughput. Nucleic Acids Res. 32, 1792–1797. https://doi.org/10.1093/nar/gkh340

Excoffier, L., Smouse, P.E., Quattro, J.M., 1992. Analysis of Molecular Variance Inferred From Metric Distances Among DNA Haplotypes: Application to Human Mitochondrial DNA Restriction Data Laurent. Genetics 131, 479–491. https://doi.org/10.3354/meps198283

Fraser, C., Casjens, S., Huang, W., 1997. Genomic sequence of a Lyme disease spirochaete, Borrelia burgdorferi. Nature 390, 580–586. https://doi.org/doi:10.1038/37551

Gatzmann, F., Metzler, D., Krebs, S., Blum, H., Sing, A., Takano, A., Kawabata, H., Fingerle, V., Margos, G., Becker, N.S., 2015. NGS population genetics analyses reveal divergent evolution of a Lyme Borreliosis agent in Europe and Asia. Ticks Tick. Borne. Dis. 6, 344–351.

Gelman, A., Su, Y.-S., 2016. arm: Data Analysis Using Regression and Multilevel/Hierarchical Models.

Haas, B.J., Delcher, A.L., Wortman, J.R., Salzberg, S.L., 2004. DAGchainer: A tool for mining segmental genome duplications and synteny. Bioinformatics 20, 3643–3646. https://doi.org/10.1093/bioinformatics/bth397

Hudson, R., Kaplan, N., 1985. Statistical properties of the number of recombination events in the history of a sample of DNA sequences. Genetics 111, 147–164.

Hudson, R., Slatkin, M., Maddison, W.P., 1992. Estimation of levels of gene flow from DNA sequence data. Genetics 132, 583–589. https://doi.org/10.1093/genetics/132.2.583

Kamvar, Z.N., Brooks, J.C., Grünwald, N.J., 2015. Novel R tools for analysis of genome-wide population genetic data with emphasis on clonality. Front. Genet. 6, 208. https://doi.org/10.3389/fgene.2015.00208

Katoh, K., Misawa, K., Kuma, K., Miyata, T., 2002. MAFFT: a novel method for rapid multiple sequence alignment based on fast Fourier transform. Nucleic Acids Res. 30, 3059–3066.

Katoh, K., Standley, D.M., 2013. MAFFT multiple sequence alignment software version 7: Improvements in performance and usability. Mol. Biol. Evol. 30, 772–780. https://doi.org/10.1093/molbev/mst010

Kurtz, S., Phillippy, A., Delcher, A.L., Smoot, M., Shumway, M., Antonescu, C., Salzberg, S.L., 2004. Versatile and open software for comparing large genomes. Genome Biol. 5, R 12.

Lemey, P., Rambaut, A., Drummond, A.J., Suchard, M.A., 2009. Bayesian phylogeography finds its roots. PLoS Comput. Biol. 5. https://doi.org/10.1371/journal.pcbi.1000520

Margos, G., Gatewood, A.G., Aanensen, D.M., Ra Hanincová, K., Terekhova, D., Vollmer, S.A., Cornet, M., Piesman, J., Donaghy, M., Bormane, A., Hurn, M.A., Feil, E.J., Fish, D., Casjens, S., Wormser, G.P., Schwartz, I., Kurtenbach, K., Hanincová, K., Terekhova, D., Vollmer, S.A., Cornet, M., Piesman, J., Donaghy, M., Bormane, A., Hurn, M.A., Feil, E.J., Fish, D., Casjens, S., Wormser, G.P., Schwartz, I., Kurtenbach, K., Ra Hanincová, K., Terekhova, D., Vollmer, S.A., Cornet, M., Piesman, J., Donaghy, M., Bormane, A., Hurn, M.A., Feil, E.J., Fish, D., Casjens, S., Wormser, G.P., Schwartz, I., Kurtenbach, K., 2008. MLST of housekeeping genes captures geographic population structure and suggests a European origin of Borrelia burgdorferi. Proc. Natl. Acad. Sci. U. S. A. 105, 8730–8735.

Nakagawa, S., Schielzeth, H., 2010. Repeatability for Gaussian and non-Gaussian data: A practical guide for biologists. Biol. Rev. 85, 935–956.

Nei, M., 1987. Molecular Evolutionary Genetics. Columbia University Press, New York.

Overbeek, R., Olson, R., Pusch, G.D., Olsen, G.J., Davis, J.J., Disz, T., Edwards, R.A., Gerdes, S., Parrello, B., Shukla, M., Vonstein, V., Wattam, A.R., Xia, F., Stevens, R., 2014. The SEED and the Rapid Annotation of microbial genomes using Subsystems Technology (RAST). Nucleic Acids Res. 42, D206–D214. https://doi.org/10.1093/nar/gkt1226

Paradis, E., 2010. Pegas: An R package for population genetics with an integrated-modular approach. Bioinformatics 26, 419–420. https://doi.org/10.1093/bioinformatics/btp696

Pfeifer, B., Wittelsbürger, U., Ramos-Onsins, S.E., Lercher, M.J., 2014. PopGenome: An efficient swiss army knife for population genomic analyses in R. Mol. Biol. Evol. 31, 1929–1936. https://doi.org/10.1093/molbev/msu136

Pollack, R.J., Telford, S.R., Spielman, A., 1993. Standardization of medium for culturing Lyme disease spirochetes. J. Clin. Microbiol. 31, 1251–1255. https://doi.org/10.1128/jcm.31.5.1251-1255.1993

Preac-Mursic, V., Wilske, B., Schierz, G., 1986. European Borrelia burgdorferi Isolated from Humans and Ticks Culture Conditions and Antibiotic Susceptibility. Zentralblatt für Bakteriol. Mikrobiol. und Hyg. 263, 112–118.

R Core Team, 2019. R: A language and environment for statistical computing. R Foundation for Statistical Computing, Vienna, Austria.

Rambaut, A., Drummond, A.J., Xie, D., Baele, G., Suchard, M.A., 2018. Posterior summarization in Bayesian phylogenetics using Tracer 1.7. Syst. Biol. 67, 901–904. https://doi.org/10.1093/sysbio/syy032

Revell, L.J., 2012. phytools: an R package for phylogenetic comparative biology (and other things). Methods Ecol. Evol. 3, 217–223. https://doi.org/https://doi.org/10.1111/j.2041-210X.2011.00169.x

Rollins, R.E., Yeyin, Z., Wyczanska, M., Alig, N., Hepner, S., Fingerle, V., Margos, G., Becker, N.S., 2021. Spatial variability in prevalence and genospecies distributions of Borrelia burgdorferi sensu lato from ixodid ticks collected in southern Germany. Ticks Tick. Borne. Dis. 12. https://doi.org/10.1016/j.ttbdis.2020.101589

Schwartz, I., Margos, G., Casjens, S.R., Qiu, W.G., Eggers, C.H., 2021. Multipartite genome of Lyme disease Borrelia: Structure, variation and prophages. Curr. Issues Mol. Biol. 42, 409–454. https://doi.org/10.21775/cimb.042.409

Tajima, F., 1989. Statistical Method for Testing the Neutral Mutation Hypothesis by DNA Polymorphism. Genetics 123, 585–595.

Takano, A., Toyomane, K., Konnai, S., Ohashi, K., Nakao, M., Ito, T., Andoh, M., Maeda, K., Watarai, M., Sato, K., Kawabata, H., 2014. Tick Surveillance for Relapsing Fever Spirochete Borrelia miyamotoi in Hokkaido, Japan. PLoS One 9, e104532.

Tavaré, S., 1986. Some Probabilistic and Statistial Problems in the Analysis of DNA Sequences. Lect. Math. teh Life Sci. 17, 57–86.

Ullrich, K., 2020. CRBHits: From Conditional Reciprocal Best Hits to Codon Alignments and Ka/Ks in R. J. Open Source Softw. 5, 2424. https://doi.org/10.21105/joss.02424

Wallace, R.G., HoDac, H., Lathrop, R.H., Fitch, W.M., 2007. A statistical phylogeography of influenza A H5N1. Proc. Natl. Acad. Sci. 104, 4473–4478. https://doi.org/10.1073/pnas.0700435104

Walter, K.S., Carpi, G., Caccone, A., Diuk-Wasser, M.A., 2017. Genomic insights into the ancient spread of Lyme disease across North America. Nat. Ecol. Evol. 1, 1569–1576. https://doi.org/10.1038/s41559-017-0282-8