DOI: 10.5061/dryad.m37pvmd28
title: Out of Asia? Expansion of Eurasian Lyme borreliosis causing genospecies display unique evolutionary trajectories

description: 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.
description: <em>Isolates, culturing, and DNA
 extraction</em> The major aim of this study was
 to estimate the direction of colonization in each
 <em>Borrelia</em> 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
 <em>Borrelia</em> isolates coming from <em>B</em>.
 <em>afzelii</em> (Asian, n=20; European, n=13),
 <em>B</em>. <em>bavariensis</em> (Asian, n =27;
 European, n=19), and <em>B</em>. <em>garinii</em>
 (Asian, n=25; European, n=32). Of these, 52 are novel
 <em>Borrelia</em> 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
 (<em>B</em>. <em>afzelii</em>, n=11;
 <em>B</em>. <em>garinii</em>, n=25; and
 <em>B</em>. <em>bavariensis</em>, n=19) and 12
 Japanese isolates (<em>B</em>.
 <em>bavariensis</em>, n=8; <em>B</em>.
 <em>garinii</em>, n=4) were provided by the German National
 Reference Center for <em>Borrelia</em> 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 <em>B</em>. <em>bavariensis</em> (n=6)
 (Becker et al., 2020) and <em>B</em>.
 <em>garinii</em> (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. <em>Borrelia</em> 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<sup>8</sup> cells per mL, at which point whole genomic DNA
 was extracted. Genomic DNA from all European isolates was extracted using
 a Maxwell<sup>®</sup> 16 LED DNA kit (Promega, Madison, WI,
 USA) and from all Japanese and Russian isolates using the
 Wizard<sup>®</sup> Genomic DNA Purification Kit (Promega, WI,
 USA). DNA quality (260/280) and concentration were measured using a
 NanoDrop<sup>®</sup> 1000 photometer (Thermo Fisher
 Scientific, Waltham, MA, USA) and a Qubit<sup>®</sup> 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. <em>Whole genome sequencing, assembly, and determining orthologus genes</em> Within the course of this study, we aimed to reconstruct the evolutionary history of <em>B</em>. <em>afzelii</em>, <em>B</em>. <em>bavariensis</em>, and <em>B</em>. <em>garinii</em> 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 <em>Borrelia</em> genome and the majority of the sequence (&gt;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 <em>Borrelia</em> genomes (Becker et al., 2020). Pacific Bioscience sequences were obtained for three <em>B</em>. <em>bavariensis</em> isolates (PBi, A104S, and NT24) (Becker et al., 2020) and three <em>B</em>. <em>garinii</em> isolates (PHeI, PBr, and NT31; see Suppl. Met.). Additionally, three <em>B</em>. <em>afzelii</em> 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 <em>B</em>. <em>bavariensis</em> 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 <em>crb2rbh</em> 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 <em>igraph</em> 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 <em>pegas</em> (Paradis, 2010) as well as the proportion of isolates carrying the synteny member group.  <em>Recombination analysis</em> Recombination is known to be low on the <em>Borrelia</em> 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 <em>x<sub>i</sub></em> is the average violation of block <em>i</em> over all comparisons, <em>μ<sub>within</sub></em> is the average within-block violation score, and <em>sd<sub>within</sub></em> 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: <em>blastn</em>) 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 <em>phytools</em> (Revell, 2012). <em>Phylogenetic reconstruction including biogeographical inference</em> In the aim to estimate the pattern of <em>Borrelia</em> 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. <em>Population genetic analyses</em> 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 (<em>π</em>) (Nei, 1987) and Tajima’s <em>D</em> test statistic (Tajima, 1989) were estimated in the package <em>pegas</em> (Paradis, 2010). Analysis of molecular variance (AMOVA) (Excoffier et al., 1992) was performed using the package <em>poppr</em> (Kamvar et al., 2015) whereas <em>F<sub>ST</sub></em> (Nei, 1987) and <em>D<sub>XY</sub></em> (Hudson, Slatkin, & Maddison, 1992) were estimated with the package PopGenome (Pfeifer et al., 2014). <em>Identification of plasmid content through plasmid partitioning genes</em> The primary aim of this study was to utilize full chromosome sequences to estimate the direction of <em>Borrelia</em> colonization between Europe and Asia as well as quantify the variability in genetic diversity between populations. Even so, the <em>Borrelia</em> genome contains many accessory plasmid types which have also been shown to be influenced by demographic events such as the bottleneck observed in <em>B</em>. <em>bavariensis</em> 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 <em>Borrelia</em> 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: <em>blastn</em>) 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 <em>t</em>.<em>test</em> from the base R package (R Core Team, 2019). Classical multidimensional scaling (MDS) was run using the <em>cmdscale</em> 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 <em>glmer</em> function from the package <em>lme4</em> (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 <em>sim</em> function from the package <em>arm</em> (Gelman and Su, 2016). Residual error was calculated according to Nakagawa & Schielzeth (2010). <em>References</em> 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 <em>Borrelia bavariensis</em>. 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 <em>Borrelia burgdorferi</em>, 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 <em>Borrelia burgdorferi</em> 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, <em>Borrelia burgdorferi</em>. 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 <em>Borrelia burgdorferi</em>. 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 <em>Borrelia burgdorferi</em> 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 <em>Borrelia burgdorferi</em> 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 <em>Borrelia</em>: 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