The hare syphilis agent is related to, but distinct from, the treponeme causing rabbit syphilis

The treponemes infecting lagomorphs include Treponema paraluisleporidarum ecovar Cuniculus (TPeC) and ecovar Lepus (TPeL), infecting rabbits and hares, respectively. In this study, we described the first complete genome sequence of TPeL, isolate V3603-13, from an infected mountain hare (Lepus timidus) in Sweden. In addition, we determined 99.0% of the genome sequence of isolate V246-08 (also from an infected mountain hare, Sweden) and 31.7% of the genome sequence of isolate Z27 A77/78 (from a European hare, Lepus europeaus, The Netherlands). The TPeL V3603-13 genome had considerable gene synteny with the TPeC Cuniculi A genome and with the human pathogen T. pallidum, which causes syphilis (ssp. pallidum, TPA), yaws (ssp. pertenue, TPE) and endemic syphilis (ssp. endemicum, TEN). Compared to the TPeC Cuniculi A genome, TPeL V3603-13 contained four insertions and 11 deletions longer than three nucleotides (ranging between 6 and2,932 nts). In addition, there were 25 additional indels, from one to three nucleotides long, altogether spanning 36 nts. The number of single nucleotide variants (SNVs) between TPeC Cuniculi A and TPeL V3603-13 were represented by 309 nucleotide differences. Major proteome coding differences between TPeL and TPeC were found in the tpr gene family, and (predicted) genes coding for outer membrane proteins, suggesting that these components are essential for host adaptation in lagomorph syphilis. The phylogeny revealed that the TPeL sample from the European brown hare was more distantly related to TPeC Cuniculi A than V3603-13 and V246-08.


Introduction
In addition to human pathogens such as Treponema pallidum ssp.pallidum (TPA), ssp.pertenue (TPE), and ssp.endemicum (TEN), causing syphilis, yaws, and bejel, respectively-the Treponema genus also contains other animal pathogens [1,2].Lagomorph syphilis is caused by Treponema paraluisleporidarum [3], with its ecovar Cuniculus (TPeC) and ecovar Lepus (TPeL), infecting rabbits and hares, respectively.Current knowledge centers almost exclusively around the TPeC strain Cuniculi A, isolated in 1939 in Maryland, USA [4], which is still maintained in the laboratory.TPeC strain Cuniculi A is genetically very related to the human pathogen T. pallidum (98% genome identity on the DNA level; [2,5]).This study uses recently introduced name Treponema paraluisleporidarum [3], as this is consistent with the literature use in the last decade.However, since the GenBank database uses only validated names, the sequences were submitted as Treponema paraluiscuniculi.(T.paraluiscuniculi was validated decades ago and would not be eligible for validation according to current rules since it is highly related to Treponema pallidum [6].Like human syphilis, TPeC causes sexually transmitted infections and known to be associated with skin ulcers in the anogenital region, the face, and/ or the paws of rabbits (Oryctolagus cuniculus f. domestica) [7].Sexually transmitted TPeL infections in European brown hares (Lepus europaeus) and mountain hares (Lepus timidus), are in most cases clinically inapparent or associated with orofacial and anogenital proliferative crusty skin lesions occurring at mucocutaneous junctions [7][8][9].
Based on serological screening studies that analyzed samples from more than 1500 animals, the average seroprevalence of syphilis in wild European lagomorphs ranges from 28% to 66%, depending on the serological test(s) used and the sample's country of origin [10][11][12][13][14].The prevalence of hare syphilis has no clear geographical gradient, however, it is known that seroprevalence against TPeL negatively correlates with altitude of sampling areas [11].
Genetic studies on the causative agent of lagomorph syphilis are limited and almost exclusively based on a single available genome sequence of the TPeC strain Cuniculi A, which is composed of circular chromosome of 1,133,390 bp with no plasmids [6].So far, a single study [3] tested the genetic relatedness between TPeL and TPeC.The study revealed four nucleotide differences in the 2002 bp-long sequence obtained from partial sequences of the 16S rRNA gene, the DNA region downstream of the 16S rRNA gene, and the sequence within the TPCCA_0225 gene including one nt difference found within the 16S rRNA genes [3].
Based on infection experiments conducted by Lumeij et al. [3] and on the genetic relatedness between TPeC and TPA [5,6], we predicted that hare-infecting strains would be phylogenetically ancestral to rabbit-infecting strains and would be highly related to TPeC.In this study, we determined the complete genome sequence of Treponema paraluisleporidarum ecovar Lepus, isolate V3603-13.We used the previously described Pooled Sequence Genome Sequencing (PSGS) technique to compare the V3603-13 genome to the once-determined genomic sequence TPeC strain Cuniculi A. We found that the treponeme causing hare syphilis was in many ways similar but distinct from the treponeme causing rabbit syphilis.

Lagomorph samples
We obtained samples from mountain hares (Lepus timidus) hunted in Sweden.TPeL sample V3603-13 was obtained from a hare sampled close to Enånger in county Ga ¨vleborg (latitude 61.32470, longitude 17.02430) in 2013.Sample V246-08 came from a hare sample hunted in 2008 near Orsa in the county of Dalarna (latitude 61.30180, longitude 14.28340).Both samples were taken from hares submitted within the Swedish general (passive) wildlife disease surveillance program, where hunters or anybody else can report findings of disease or mortality in wildlife.The hares were shot by hunters and when they noted scabs and lesions on the lips and nose, they contacted Swedish Veterinary Agency (SVA).Upon request from SVA, the carcasses were sent to the institute for diagnostic examinations.No further information is available on the location of lesions.The brown hare isolate, Z27 A77/78, was collected in The Netherlands in 2010 and was described in an earlier study [3].

Sampling and extraction of DNA
Samples from the skin lesions (full skin thickness tissue samples) were taken from dead animals and stored in a -80 degree Celsius biobank freezer until thawed and used in this study.DNA was extracted from tissue material using the QIAamp DNA Mini Kit (Qiagen, Hilden, Germany) following the manufacturer's protocol with some minor modifications.Briefly, we extracted DNA from tissue samples according to the protocols published by Hisgen et al. [11].Subsequently, glycogen precipitation was performed to clean and concentrate the DNA as described by Knauf et al. [15].We measured the DNA yield using a NanoDrop photometer (Thermo Fisher Scientific, Darmstadt, Germany).

Amplification of genomic DNA
The genomic DNA of all three samples was amplified using the multiple displacement amplification approach (REPLI-g kit, QIAGEN, Valencia, CA, USA).Samples were diluted one to 50 times (depending on the positivity of treponemal-specific PCR amplification from undiluted and diluted samples) and then used as a template for PSGS, as previously described [16][17][18].The DNA was amplified using PrimeSTAR GXL DNA Polymerase (Takara Bio Inc., Otsu, Japan) with 133 specific primer pairs (published by Mikalova ´et al. [19], with minor modifications) to obtain overlapping PCR products covering the entire chromosome.PCR products were amplified using touchdown PCR under the following cycling conditions: initial denaturation at 94˚C for 1 min; 8 cycles: 98˚C for 10 s, 68˚C for 15 s (the annealing temperature was gradually reduced at a rate of 1˚C per cycle), and 68˚C for 6 min; 35 cycles: 98˚C for 10 s, 61˚C for 15 s, and 68˚C for 6 min (43 cycles in total); followed by the final extension at 68˚C for 7 min.PCR products were subsequently purified using QIAquick PCR Purification Kits (QIA-GEN, Valencia, CA, USA) and mixed in equimolar amounts to produce four distinct pools, which were then used for whole genome sequencing.

Whole genome sequencing and assembly of the genomes
The sequencing library for TPeL strain V3603-13 was prepared from 1 ng of RNA-free genomic DNA using a Nextera XT DNA Sample Preparation kit (Illumina, CA, USA).Sequencing was performed using the MiSeq Reagent Kit v3 in the MiSeq system (Illumina, USA) at the Veterinary Research Institute sequencing facility (Brno, Czech Republic).The sequencing reads of individual pools were handled separately and assembled de novo with SeqMan NGen v4.1.0software (DNASTAR, Madison, WI, USA) using the default parameters.Contigs obtained from TPeL strains were then aligned to the genome of the TPeC strain Cuniculi A (CP002103.1) using Lasergene software (DNASTAR, Madison, WI, USA).Missing parts of genomes were subsequently Sanger sequenced (GATC Biotech, Germany).In addition, Min-ION sequencing (Oxford Nanopore Technologies (ONT), Oxford, UK) was used to sequence paralogous gene regions, mainly those containing tpr genes.
Both the V246-08 and Z27 A77/78 genomes were amplified using the PSGS method and then sequenced using the Illumina platform as described for the V3603-13 sample.The quality of the raw reads was checked using FastQC [20].Raw reads were trimmed in length (a minimal length of 35 nt) and quality (Phred quality score � 20) using Cutadapt [21].Preprocessed reads were mapped using the BWA MEM algorithm [22] onto the TPeL Cuniculi A genome (GenBank Acc.No. CP002103.1).The Cuniculi A reference genome was divided into four parts representing the "pools" described above.The post-processing of the reads' mapping was performed using Samtools [22], Picard (https://broadinstitute.github.io/picard/),GATK [23], and NGSUtils/bamutils [24].Low-quality mappings were omitted from analyses (i.e., mapping quality; MAPQ < 40) as well as reads mapping to multiple sites (e.g., in repetitive and paralogous regions).A minimal alignment length was set to 35 bps, the maximum number of allowed mismatches was set to 5 (or 5% of the read length), and the maximum soft-clipping was set to 5% of the read length.Consensus sequences were determined using variants detected using Samtools "variant calling" followed by Vcfutils [22] and Seqtk (https://github.com/lh3/seqtk).The breadth and depth of coverages were calculated using GATK software.

Construction of phylogenetic tree
To determine evolutionary relationships among pathogenic treponemes, corresponding parts of the available treponemal genomes (orthologous sequences) were used and included TPeC Cuniculus strain Cuniculi A (CP002103.1),TPeC Cz-2020 (MW323408), and TPA strain SS14 (CP004011.1).The phylogenetic tree was constructed using the Maximum Likelihood method based on the Tamura-Nei model and MEGA 7 software, which was also used for model testing [25].Bootstrap values were equal to 1000.

Nucleotide sequence accession number
The complete genome sequence of the TPeL V3603-13 isolate was deposited in GenBank under accession number CP097901, Bioproject PRJNA606433.Since Treponema paraluisleporidarum ecovar Lepus is not validated, the sequences were submitted as Treponema paraluiscuniculi into the GenBank database.This is the reason why sequence CP097901 was submitted as Treponema paraluiscuniculi strain L2, which was the original designation of the isolate.The sequencing reads assembled into draft genome sequences of the TPeL V246-08 (JAMZQX000000000) and Z27 A77/78 isolates were deposited in GenBank under Bioproject accession numbers PRJNA837879 and PRJNA837886, respectively.

Sequencing of TPeL isolates
TPeL isolates were sequenced using PSGS, as described earlier [16,17].Briefly, DNA was amplified with 133 pairs of specific primers to obtain overlapping PCR products [18]).To facilitate sequencing of paralogous genes containing repetitive sequences, PCR products were mixed in equimolar amounts into four distinct pools.The PCR products constituting each pool were labeled with multiplex identifier (MID) adapters and Illumina sequenced.Pools

The overall genome structure of TPeL isolate V3603-13
The complete genome sequence of TPeL isolate V3603-13 comprised 1,132,489 nucleotides with no sequencing gaps or ambiguities.The genome sequence also includes variable sites in the tprK gene (TP0897), where the most common variant (with a frequency greater than 50%) was used for the final whole genome consensus sequence.The TPeL V3603-13 genome is 901 nt smaller than the reference genome of TPeC strain Cuniculi A (GenBank acc.number CP002103).There is overall gene synteny with the TPeC strain Cuniculi A genome and the human pathogen T. pallidum (subsp.TPA, TPE, and TEN).The TPeL V3603-13 genome showed a different rrn spacer pattern compared to the TPeC Cuniculi A genome, i.e., tRNA gene for Ala/Ile in V3603-13 compared to tRNA gene for Ile/Ala in the rabbit infecting strain of Cuniculi A. Sequences of the 5S, 16S, and 23S rRNA genes were identical in both operons and were without macrolide resistance coding single nucleotide changes A2058G or A2059G in the 23S rRNA genes [26,27].The number of 60 bp-long repetitions within the arp gene (TP0433) in TPeL V3603-13 and TPeC strain Cuniculi A genomes differed and was 19 and 21 repetitions, respectively.The number of 24 bp-long repetitions was 19 in the TP0470 in TPeL V3603-13, and 6 in TPeC strain Cuniculi A. Table 2 summarizes the overall genome parameters of TPeL isolate V3603-13.

Major sequence differences between TPeL V3603-13 and TPeC Cuniculi A
Despite the overall genome synteny of TPeL V3603-13 and TPeC strain Cuniculi A, several genetic differences differentiate the genomes, including indels and single nucleotide variants (SNVs).An overview of the genetic differences between the rabbit infecting strain, Cuniculi A, and the V3603-13 genome is shown in Fig 1 and S1 Table .While most of the indels were found within genes, a minority was found in the intergenic regions (IGR).Moreover, most of the indels in the genes did not change the reading frame leading to indels of one to few amino acids.Compared to the TPeC strain Cuniculi A genome, TPeL V3603-13 contains four insertions greater than three nucleotides: a 1,874 nt-long insertion containing a TP0126c-TP0129 fragment similar to the orthologous TPA strain Philadelphia 1 genomic sequence [28], a 79 nt insertion in the IGR between TP0545-TP0546 similar to the TPA strain Philadelphia 1 genomic sequence, a 24 nt insertion in TP0577 that is similar to the TEN strain 11q/j genomic sequence [29], and a 12 nt insertion (with undetected sequence similarity) in TP0897 (tprK) (S1 Table ).
Compared to the TPeC strain Cuniculi A genome, the TPeL V3603-13 genome contained eleven deletions greater than three nucleotides: a 2,932 nt deletion similar (but larger) to the deletion present in TEN strains [29] comprising the tprFG and tmpC (TP0319) genes, a six ntlong deletion in TP0126b, a 45 nt deletion in TP0136, a 13 nt deletion in TP0146, a six nt deletion in TP0515, a 36 nt deletion in TP0548, a six nt deletion in the IGR between TP0650-TP0652, a six nt deletion in TP0733, a 17 nt deletion in the IGR between TP0895-TP0897, a six nt deletion in TP0897 (tprK), and a six nt deletion in TP0966 (S1 Table ).In addition, there were 25 additional indels with lengths of 1-3 nucleotides, altogether spanning 36 nt (S1 Table ).Three of these changes lead to in-frame insertions or deletions.Other indels located in genes lead to frameshifts, which are shown in S1 Table .A subset of insertion/deletion differences described in S1 Table was also found in TPeL isolate V246-08 supporting lack of sequencing artifacts.
Several sequentially diverse chromosomal regions (defined as multiple SNVs) were identified in TPeL V3603-13 genome: a 59 nt long region in TP0126b showing similarity to tprK of TPA SS14, a 75 nt long region in TP0126b showing similarity to tprK of TPA X-4, a 328 nt long region in TP0136 with similarity to TP0134 of TPA Philadelphia 1, a 54 nt long region in position of TP0326 similar to bamA of TPA CZ_177zB (TPA Nichols-like strain) [30], a 36 nt long region not similar to any known treponemal sequences in position of TP0548, a 26 nt long region similar to IGR TP0650-TP0652 of TPA Philadelphia 1 in position of IGR TP0650-652, a 24 nt long region not similar to any known treponemal sequences in position of TP0733, four regions upstream and in the tprK gene (3nt, 45 nt, 41nt, 136 nt) similar to TPeC or TPA SS14 sequences, and a 15 nt long region similar to TPE Fribourg-Blanc in position of TP0966 (S1 Table ).
Moreover, frameshift mutations resulting in two new pseudogenes (TP0308a, TP1035a) and major sequence changes, including protein elongation or shortening (TP0040, TP0179, TP0279, TP0471, TP0778, TP0801) were described.In addition, another frameshift mutation in existing pseudogenes was found in two pseudogenes (TP0103, TP0146) (S1 Table ).Lastly, the mglB gene (TP0545) of TPeL V3603-13 appears to be fully functional and similar to TPE strains (due to the insertion of 79 nts in IGR between TP0545-TP0546).Similarly, a nucleotide insertion in V3603-13 TP0617a results in a gene version identical to the paralogous TP0315 gene in the V3603-13 genome.Genes TP0651 and TP0896 in V3603-13 are longer but similar to the versions present in TPE strains.
As shown in Fig 2, repeat motifs within the arp gene (TP0433) and the arrangement of the repeat sequences differ between TPeL V3603-13 and TPeC Cuniculi A genomes.

Single nucleotide sequence differences between the TPeL V3603-13 and TPeC strain Cuniculi A genomes
In addition to the described differences, 309 additional nucleotide variants are dispersed throughout the genome in TPeL V3603-13 compared to the TPeC Cuniculi A reference genome.These nucleotide differences were primarily represented by single nucleotide variants (SNVs) and less frequently by double nucleotide differences (n = 8, covering 16 nucleotides).Analysis of the SNVs between TPeL samples V3603-13 and V246-08 revealed altogether 35 single nucleotide variants indicating that majority of genetic differences between hare and rabbit syphilis is shared by both V3603-13 and V246-08 samples.

The phylogenetic relationship of lagomorph and human pathogenic Treponema
To determine the genetic relatedness of TPeL/TPeC samples, a tree was built on partial TP0548 sequences available for two TPeC strains/isolates (Cuniculi A; Cz-2020) and three TPeL isolates TPeL V3603-13, V246-08, and Z27 A77/78 (Fig 3A).Sample Cz-2020 is the only other known member of TPeC and only sequences of few loci are available.Therefore, locus TP0548 was selected as the only locus with available sequences in all 5 T. paraluisleporidarum samples.TPeL isolates clustered together but separate from TPeC isolates.To determine the genetic relatedness of both ecovars using a more robust analysis, a tree was built on partial genome sequences available for the complete genome of TPeL V3603-13 and TPeL isolates V246-08 and TPeL Z27 A77/78.Altogether, about one-third of the genome length (355,133 nucleotides) was used to construct the phylogenetic tree (Fig 3B).Both TPeL isolates V3603-13 and V246-08, which were isolated from L. timidus, clustered closely with the TPeC Cuniculi A genome.At the same time, the Z27 A77/78 sample, taken from the L. europeus in The Netherlands [3], was more distantly related.

Discussion
In this study, we have determined the first whole genome sequence of the causative agent of hare syphilis.TPeL isolate V3603-13 was obtained from an infected Swedish mountain hare.In addition, we generated two partial genome sequences for TPeL isolate V246-08, also from an infected mountain hare in Sweden, and Z27 A77/78 that came from a European brown hare in The Netherlands [3].While the sample quality of V246-08 and Z27 A77/78 (quality refers to the amount of available treponemal DNA) did not allow us to sequence the whole genome, considerable parts of the corresponding genomes were sequenced with the shortest  [32,33].There was a total of 777 positions in the final dataset.B. The tree was built on available genome sequences.There were 355,133 positions in the final dataset.The trees were constructed using the Maximum Likelihood method based on the Tamura-Nei model.Bootstrap support is shown next to the branches.The scale shows the number of substitutions per site.As an outgroup, the TPA SS14 sequence was used.All positions containing gaps and missing data were omitted.https://doi.org/10.1371/journal.pone.0307196.g003determined sequence from Z27 A77/78 being equal to about a third of the TPeC strain Cuniculi A genome, i.e., 31.7%,355,133 nucleotides.
Until now, only a limited number of studies have been published on the genetics of TPeL [3,34] even though infections caused by TPeL appear to be widespread among European hares [11][12][13][14]35].Since most of the previous and current TPeL infections are subclinical and detectable only with serology and/or direct detection of TPeL DNA [8], the scarcity of clinical symptoms could explain the lack of scientific data on treponemal infections in hares.
The genome of TPeL V3603-13 comprised of 1,132,489 nt and is the smallest genome determined to date of the group of closely related pathogenic treponemes, including human pathogenic TPA, TPE, and TEN.At the whole genome level, genetic similarity between TPeL and TPeC was > 99.8%, with almost complete gene synteny between the genomes.Interestingly, both the TPeL V3603-13 and TPeC Cuniculi A genomes showed different rrn spacer patterns (Ala/Ile and Ile/Ala, respectively), a feature previously noticed within TPA and TPE strains as a result of reciprocal intragenomic translocation events [2,34].The number of repetitions within the arp gene (TP0433) in TPeL V3603-13 was less than TPeC Cuniculi A (19 vs. 21), and the gene was composed of different sequence motifs compared to the arp of Cuniculi A (Fig 2) [31].Among TPA isolates, the number and the structure of the arp repeats were shown to generally correlate with whole genome phylogeny [36].Interestingly, the Cuniculi A strain arp gene analyzed by Strouhal et al. [5] and S ˇmajs et al. [6] differed in the sequence of repeat motifs from the same strain analyzed by Harper et al. [31], suggesting that these gene repeat components are prone to frequent reshuffling and accumulation of point mutations.Consistent with this observation, previous studies on the number of arp repeat motifs among TPA isolates revealed differences in the repeat numbers between whole blood and swab samples taken from the same patient [37], again suggesting the genetic plasticity of this locus.Both the TPeC and TPeL arp genes contained multiple repeat motifs (Cuniculi A, n = 4; V3603-13, n = 5), a feature that has only been described for TPA and TPeC strains but not for TPE (including stains infecting nonhuman primates [38] and TEN strains [31]. One major difference between the rabbit infecting TPeC Cuniculi A and the hare TPeL V3603-13 genome was the presence of an 1874 nt-long insertion in V3603-13.This insertion mainly corresponded to the TP0126c-TP0129 genes in other TPA and TPE treponemes, which suggests that the corresponding region was deleted during the evolution of the Cuniculi A strain.In a subpopulation of the TPA Nichols strain, a 1.3-kb deletion at this genomic locus (i.e., in TP0126) was found in an earlier study [39]; both the TPA Nichols and TPeC Cuniculi A deletions overlapped by 376 nt.Another previous study revealed that TP0126, TP0126b, TP0126c, TP0127d, TP0128, and TP0130 (in addition to other genes) have a modular genetic structure enabling rapid genetic diversification of treponemal strains [40].At the same time, these regions, due to the presence of repetitive sequences, may also have less genetic stability and be prone to deletion events.Since the Nichols strain has been propagated in rabbits for some time, the deletion of this region may represent one of the adaptations to infecting rabbits (or hares).At least two additional insertions (in TP0545-546 and TP0577) in the V3603-13 genome matched sequences present in other TPA and TPE strains, suggesting additional deletions in the evolution of the rabbit pathogen, TPeC.Based on cross-infection experiments [3], it was predicted that TPeL can infect both rabbits and hares while TPeC infects only rabbits.Therefore, the deleted regions in TPeC may be important relative to hare infections.However, genome manipulation experiments would be needed to prove this and better understand the host-pathogen evolution of treponemal pathogens.
Another striking feature of the TPeL V3603-13 genome is a 2932 nt-long deletion in V3603-13 affecting tprFG and tmpC (TP0319) genes (compared to TPeC strain Cuniculi A).The genome of TPeC Cuniculi A has a similar but smaller deletion in this region (compared to TPA strain Nichols), suggesting the independent origin of these deletions.Moreover, the surrounding genes, including TP0309-311, TP0313, TP0315, and TP0318, contain frameshift mutations in both TPeC Cuniculi A and TPeL V3603-13 genomes, and it is therefore possible that the entire region is not needed for infection of lagomorphs, including both rabbits and hares.Yet, it is not clear if the tprFG region is required for infection of humans since similar but smaller deletions in the TEN strains were found in subpopulations of individual TEN strains that were experimentally propagated in rabbits [17,41].On the other hand, the TEN strains isolated directly from Cuban patients suspected of having syphilis [42] contained a deletion at the tprFG locus [43] which was similar to one described in the TEN reference strains Bosnia A and Iraq B [17,41].This suggests that the deletion does not preclude infection in humans nor the emergence of early syphilis-like symptoms [29,42,44].
Analysis of the TP0136 gene in the V3603-13 genome revealed a 328 nt-long region similar to TP0134 in TPA Philadelphia 1 [28] while the TP0136 gene in Cuniculi A had a sequence similar to TP0133 (Fig 4).Since the TP0134 gene (including the donor sequence) was deleted from both the TPeC Cuniculi A and TPeL V3603-13 genomes, this finding cannot be explained by sequence gene transfer from the TP0134 locus.An explanation is the presence of a common ancestor, i.e., for the TPeC Cuniculi A and TPeL V3603-13 strain, with a TP0134like sequence at the TP0136 locus (Fig 4).Subsequently, a deletion might have occurred.The lineage then became ancestral to the TPeL genome while additional gene conversion leading to replacement of the TP0134-like sequence by the TP0133-like sequence occurred during the evolution of the rabbit infecting TPeC genome of Cuniculi A (Fig 4).Although this model suggests TPeC being an evolutionarily modern version of TPeL, other alternative explanations, including gene recombination events [45] are also possible.
The TP0136 locus encodes for the fibronectin-binding protein [46] and is known to be a recombinant and positively selected locus in TPA treponemes [28,47].Moreover, this locus has been predicted to have a modular structure that can lead to relatively frequent sequence changes and indels [40,48].However, the part of TP0136 with the modular structure was not affected during the above-described replacements (Fig 4).All this evidence points to the importance of the TP0136 locus during human and animal infections.Besides deletions, insertions, and stretches of nucleotide sequences resulting from gene conversions, we found 309 single nucleotide variants (SNVs) between the TPeC Cuniculi A and TpeL V3603-13 genomes.These differences were present in single or two nt sequences (n = 8), suggesting that these events represented accumulated mutations.Although mutation rates in TPeL and TPeC remain unknown, previous studies estimated the mutation rates of the genetically related treponemes TPA and TPE.The mutation rate was estimated to be 2.8-4.1 × 10 −10 per site per generation [49,50], which corresponds to a mutation rate of 0.846-1.21× 10 −7 per nucleotide site per year.Other studies, based on temporal analyses with relaxed clock models, estimated slightly higher treponemal mutation rates, i.e., 3.02 × 10 −7 [51] and 6.6 × 10 −7 [52] per nucleotide site per year, which would correspond to 0.5-3 kiloyears of separate evolution for lagomorph treponemes.We note here that these mutation rate estimates, based on TPA/E strains in humans, do not necessarily correspond to lagomorph infecting TPeC/L.Moreover, mutation rates could vary over time, especially during the early phases of adaptation to a new host (species), and mutation rates could be even lower in TPA/TPE [18,50].
The mutations causing resistance to macrolide antibiotics (i.e., A2058G or A2059G mutations in the 23S rRNA genes [26,27,53]) were not found among the analyzed TPeL isolates.TPeL strains and TPeC strain Cuniculi A remain the only group of pathogenic treponemes where these mutations have not been detected, in contrast to TPA [26,27], TPE [54,55], and TEN [56,57].The absence of macrolide resistance in TPeL and TPeC is consistent with the absence of macrolide resistance in TPE isolated from wild nonhuman primates [15,38,58,59].Although speculative and theoretical resistance could quickly emerge in wildlife, too, particularly under selection pressure caused by the mass administration of antimicrobial treatments.

Conclusions
We determined the first complete genome sequence of Treponema paraluisleporidarum ecovar Lepus (TPeL V3603-13) isolated from a naturally infected mountain hare in Sweden and provided draft genome sequences of two additional TPeL strains.The agent causing hare syphilis was found to be similar but distinct from the rabbit syphilis treponeme; the major predicted proteome differences were associated with the Tpr proteins and outer membrane proteins/ antigens.Based on previous estimations of TPA and TPE mutation rates, both TPeL and TPeC appear to be separated by 0.5-3 kiloyears of lagomorph treponeme evolution.However, more samples from naturally infected hares and rabbits need to be analyzed to understand better the genetic diversity of lagomorph syphilis and the phylogeny of lagomorph treponemes.

Fig 1 .
Fig 1.An overview of genome differences between TPeL V3603-13 and TPeC Cuniculi A. Insertions and deletions are shown above and below the schematic representation of the TPeL isolate V3603-13 chromosome, respectively.Green areas represent sequentially diverse regions, and the lengths of these regions are shown in green letters, while indel lengths are shown in black letters.https://doi.org/10.1371/journal.pone.0307196.g001

Fig 3 .
Fig 3. Molecular phylogeny of TPeL and TPeC samples.A. The tree was built on partial TP0548 sequences (TPeC Cz-2020 and TPA SS14)[32,33].There was a total of 777 positions in the final dataset.B. The tree was built on available genome sequences.There were 355,133 positions in the final dataset.The trees were constructed using the Maximum Likelihood method based on the Tamura-Nei model.Bootstrap support is shown next to the branches.The scale shows the number of substitutions per site.As an outgroup, the TPA SS14 sequence was used.All positions containing gaps and missing data were omitted.