Molecular Characterization of Borrelia persica, the Agent of Tick Borne Relapsing Fever in Israel and the Palestinian Authority

The identification of the Tick Borne Relapsing Fever (TBRF) agent in Israel and the Palestinian Authority relies on the morphology and the association of Borrelia persica with its vector Ornithodoros tholozani. Molecular based data on B. persica are very scarce as the organism is still non-cultivable. In this study, we were able to sequence three complete 16S rRNA genes, 12 partial flaB genes, 18 partial glpQ genes, 16 rrs-ileT intergenic spacers (IGS) from nine ticks and ten human blood samples originating from the West Bank and Israel. In one sample we sequenced 7231 contiguous base pairs that covered completely the region from the 5′end of the 16S rRNA gene to the 5′end of the 23S rRNA gene comprising the whole 16S rRNA (rrs), and the following genes: Ala tRNA (alaT), Ile tRNA (ileT), adenylosuccinate lyase (purB), adenylosuccinate synthetase (purA), methylpurine-DNA glycosylase (mag), hypoxanthine-guanine phosphoribosyltransferase (hpt), an hydrolase (HAD superfamily) and a 135 bp 5′ fragment of the 23S rRNA (rrlA) genes. Phylogenic sequence analysis defined all the Borrelia isolates from O. tholozani and from human TBRF cases in Israel and the West Bank as B. persica that clustered between the African and the New World TBRF species. Gene organization of the intergenic spacer between the 16S rRNA and the 23S rRNA was similar to that of other TBRF Borrelia species and different from the Lyme disease Borrelia species. Variants of B. persica were found among the different genes of the different isolates even in the same sampling area.


Introduction
Tick Borne Relapsing Fever (TBRF) is characterized by recurring fever attacks, usually of decreasing intensity and accompanied by nonspecific symptoms like myalgia, headache and gastrointestinal symptoms. It is caused by a dozen different spirochetes species of the genus Borrelia that are endemic in different geographical areas [1]. Each Borrelia species is transmitted within a geographical range by a specific species of soft ticks of the genus Ornithodoros (Argasidae) [2]. In Israel and in the Palestinian Authority the main agent of TBRF is described in the literature as Borrelia persica transmitted by the cave tick Ornithodoros tholozani [3,4,5]. However, because B. persica has not been cultured yet ''in vitro'', its characterization is mainly based on the geographic distribution of TBRF cases, the clinical and epidemiological setting, its morphology and finally on the presence of the vector O. tholozani as a source of the transmission. Accordingly the distribution of B. persica covers Central Asia (Iran, Kazakhstan, Kyrgyzstan, Tajikistan, Turkmenistan, Uzbekistan, Afghanistan and India) [3] and in the Middle East (Iraq, Jordan, Syria, Israel, and Egypt) [3].
TBRF is a rare disease in Israel but probably under-reported partly due to the fact that the disease is relatively mild and usually without a serious outcome [3]. The diagnosis and presence of TBRF in the Palestinian Authority is even less studied although geographical conditions suggest that the disease and the bacterium-host pair (B. persica-O. tholozani) are also present. Indeed, cases recorded in Jordan between 1959 and 1969 included cases from the West Bank of the River Jordan [6].
Molecular data on B. persica are very scarce and limited to a single complete sequence of the 16S ribosomal RNA (rRNA) gene (rrs) of a Persian strain [7] and partial sequences of the 16S rRNA (rrs), of the flagellin (flaB) [2] and the glycerophosphodiester phosphodiesterase (glpQ) genes [4].
A thorough molecular study is the prerequisite for the development of molecular assays for characterization of the Borrelia species involved in relapsing fever in Israel and the Palestinian Authority and consequently for the development of specific epidemiological and diagnostic assays. In this study, we sequenced three complete 16S rRNA genes, 12 partial flaB genes, 18 partial glpQ genes, 16 rrs-ileT intergenic spacers (IGS) from nine ticks and ten human blood samples originating from the West Bank and Israel. In one sample we sequenced 7231 contiguous base pairs that covered completely the region from the 59 end of the 16S rRNA gene (rrs) to the 59end of the 23S rRNA gene (rrlA). Using this information we were able to characterize by molecular methods the etiological agent of TBRF in Israel and the Palestinian Authorities present in O. tholozani and in human samples and demonstrate that it is very close if not identical to the reference stain of B. persica [7] thus providing definite evidence that B. persica is the main agent of TBRF in the Palestinian Authority and in Israel.

Ethics Statement
The study was considered as part of a routine program for TBRF diagnosis that could improve this diagnosis for future cases of TBFR and its control. The Helsinki committee of Shaare Tzedek ruled that in this particular case, formal IRB approval and written consent from patients are not required given that medical care would not be modified by the research process, given the retrospective nature of the research and that the research did not involve any procedures for which written consent is normally required outside of the research context and because the study used samples that were routinely collected for use in approved routine tests (Waiver p 27/10). All samples were anonymously obtained, no human experimentation was conducted and no human genetic study was performed.

Samples
The TBRF diagnosis in patients was established as previously reported [8]. Samples of human blood were sent to the Parasitology Reference Center (Ministry of Health, Jerusalem) and to Israeli hospitals for TBRF diagnosis. Fresh human blood samples were examined by dark field microscopy for viable Borrelia and frozen at 280uC for molecular studies.
Collection of ticks in the West Bank was conducted by Al Quds University with the collaboration of the Palestinian Ministry of Health using CO 2 traps as described previously [2]. Specimens of ticks were collected from different caves in the West Bank and were preserved in 70% alcohol. Ticks were phenotypically identified as O. tholozani by the Entomology Laboratory Israeli Ministry of Health, Jerusalem or by the Palestinian Ministry of Health, Ramallah. The list and location of the samples investigated are given in Table 1.

DNA extraction and screening for the presence of Borrelia species
Total DNA was extracted from individual ticks or from blood samples obtained from infected patients as described by Assous et al. [2]. DNeasy blood & tissue purification kit (Qiagen, Hilden, Germany) was used for DNA extractions as recommended by the manufacturer.
Preliminary screening for the presence of Borrelia in the samples was performed by amplifying a 750 bp partial fragment of the flagellin gene (flaB) by a polymerase chain reaction (PCR) with primers BOR1 and BOR2 as described by Assous et al. [2]. Details of all primers used in this work are listed in Table 2. Two negative controls, water as well as DNA extracted from uninfected ticks were included in each run. Amplified DNA was purified using the Qiaquick PCR purification kit (Qiagen, Hilden Germany) or ExoSAP-IT (USB, Cleveland, USA) as recommended by the manufacturers and used for sequencing (Hylabs, Rehovot, Israel). The amplified fragments were sequenced directly with primers used for the amplification reaction or after cloning in the EcoRV site of plasmid pBluescript using the T7 and T3 universal primers. In all cases both strands of each fragment were sequenced. Sequences were analyzed with the Vector NTI advance 11 software (Invitrogen, UK).
Amplification and sequencing of the complete 16S rRNA (rrs) gene Overlapping fragments of the 16S ribosomal RNA gene (rrs) were amplified by PCR from the blood of patients using the following pairs of primers: 16S59-16S4R, 16S1F-16S7R, 16S2F-16S10R, 16S6F-16S11R, and 16S8F-16S39 (Table 2). One unit of Phusion DNA polymerase (Finnzymes, Keilaranta, Finland) was used for DNA amplification in a 50 ml reaction mixture containing buffer GC (supplied by the manufacturer). The PCR conditions were: initial denaturation at 98uC for 30 sec, followed by 30 amplification cycles (98uC for 10 sec, 56uC for 20 sec, 72uC for 30 sec), and a final extension step at 72uC for 6 min. A 59 end fragment (16S59-16S3R) and a 39end fragment (16S9F-16S39) were cloned in the EcoRV site of the plasmid pBluescript and sequenced using the T7 and T3 universal primers.

Amplification and sequencing of the flagellin (flaB) gene
Determination of the flagellin type of different isolates was performed by sequencing a flaB fragment (750 bp) amplified using the genus-specific set of primers BOR1 and BOR2 (Table 2). For PCR amplification 1U BIOTAQ TM DNA polymerase (Bioline GMBH Germany) was used in a reaction mixture containing 1.6 mM MgCl 2 . Reaction conditions were as follows: initial denaturation step of 95uC for 1 min, followed by 40 cycles of 94uC for 30 sec, 57uC for 30 sec and 72uC for 30 sec with a final extension step of 72uC for 5 min.

Amplification and sequencing of the glycerophosphodiester phosphodiesterase (glpQ) gene
After alignment of the glpQ gene sequences available in GenBank (see below), areas of high homology were identified and used to design primers glpqF and glpqR for the amplification and cloning of the glpQ gene from B. persica (Table 2). These primers delineate an 800 bp fragment covering 80% of the gene excluding the C-and Ntermini. The PCR conditions were as described for the rrs gene. In the case of several tick isolates the glpQ fragment was not detected or was at the limit of detection after the first round of amplification; in these cases a nested protocol was applied using the same cycling conditions with an internal set of primers: glpqF and glpqRn ( Table 2). The resulting fragments were sequenced and compared by phylogenetic analysis to glpQ genes sequences of B. crocidurae (AF247151), B. recurrentis (AF247155), B. hermsii (U40762), B. parkeri (AF247156), B. turicatae (AF247157) and B. coriaceae (AF247158).  Sequencing of the rrs-rrlA spacer region Sequencing the rrs-ileT intergenic spacer (IGS). The rrs-ileT IGS portion of the rrs-rrlA spacer region was amplified using nested PCR with primers described by Bunikis et al. [9]. The first round of amplification was performed using Phusion DNA polymerase (Finnzymes, Keilaranta, Finland) and primers IGSaF and IGSaR (Table1). The following PCR conditions were used for these reactions: initial denaturation at 98uC for 30 sec followed by 35 cycles of 98uC for 10 sec, 56uC for 30 sec, 74uC for 30 sec and a final extension step of 72uC 10 min. This was followed by a nested reaction performed using the primers IGSaFn and IGSaRn ( Table 2) and PCR conditions as described for flaB.
Sequencing of the ileT-rrlA region. By aligning sequences of intergenic regions of Borrelia species available in the GeneBank we were able to identify regions of homology that were used to design primers ( Table 2) for amplification and sequencing of four overlapping fragments covering the entire region between the Ile tRNA and the 23S rRNA genes ( Table 2). The fragment between ileT and purA was amplified by a nested PCR reaction. Primers IGSb1F and IGSb1R were used in the first round and primers IGSb1Fn and IGSb1Rn were used in the nested reaction. Two overlapping fragments, a 1.3 kb fragment amplified with primers IGSb2F and IGSb2R and a 1.2 kb fragment amplified by primers IGSb3F and IGSb3R comprised the region between the purA and the hydrolase genes. The fragment from the hydrolase to the rrlA gene was amplified by a nested reaction with primers IGSb4F and IGSb4R in the first round and IGSb4Fn and IGS4Rn for the nested reaction. PCR conditions were as described for the rrs gene. Sequencing was performed using the primers used for amplification. The sequence was completed by primer walking (Table 2).
Phylogenic analysis of the different loci was performed by using the MEGA software (version 4.1) [10] after multiple alignments of sequences by CLUSTAL W (1.83) [11]. Distance options were computed according to Maximum Composite Likelihood method [12]. Phylogenic trees were generated by Unweighted Pair-Group Method with Arithmetic averages (UPGMA) assuming that TBRF Borrelia species are ultrametric [13]. A bootstrap value of 250 replicates was taken to obtain more confidence in drawing parameters.
For the multilocus sequence analysis (MLSA), sequences from different loci (rrs-ileT IGS, flaB, glpQ, purA) of seven B. persica isolates described in this work (six from human hosts: H1015, H1201, HL2610, HS3011, H1370, H1374 and one from tick: TGd1) were compared with homologue genes from other available TBRF Borrelia species. For each gene (purA, glpQ, flaB and IGS), the orthologous nucleotide sequences from four Borrelia species (B. recurrentis, B. duttonii, B. hermsii, B. turicatae) and the B. persica isolates were aligned with CLUSTAL W (1.83) [11]. The alignments were then concatenated and the resulting long alignment was used to construct a phylogenetic tree. The tree was generated with the program PHYML (v2.4.5) [14]. Default values were used except for 100 bootstraps. The resulting tree was drawn using the Interactive Tree of Life web interface (http://itol.embl.de) [15].

Results and Discussion
The failure to cultivate in vitro B. persica has prevented a comprehensive taxonomical study of this spirochete. When starting this work, available B. persica sequences were limited to one sequence of the 16S rRNA gene (rrs) [7] and fragments of the flaB [2] and glpQ genes [4]. To study the phylogenic relationship between Israeli and Palestinian TBRF Borrelia isolates and the B. persica prototype Iranian strain [7], we sequenced several loci from isolates collected in a variety of geographical location in the West Bank and in Israel from infected humans and from ticks.
Analysis of the genomic region containing the 16S rRNA gene (rrs) and the entire intergenic spacer between the 16S rRNA gene and the 23S rRNA gene (rrs-rrlA) By identifying areas of homology among various Borrelia species, in combination with primer walking, we sequenced a 7231 bp genomic region comprising the rrs gene and the entire rrs-rrlA intergenic spacer from a human blood sample positive for Borrelia by microscopy (isolate HL2610, accession number: HM131216). This region comprised in the following order: the 1526 bp of the rrs gene and a 5705 bp region from the 16S rRNA to the 59end of the 23S rRNA genes, containing the sequences encoding the Ala tRNA (alaT), the Ile tRNA (ileT), the adenylosuccinate lyase (purB), the adenylosuccinate synthetase (purA,), the methylpurine-DNA glycosylase (mag), the hypoxanthine-guanine phosphoribosyltransferase (hpt), the hydrolase (HAD superfamily) genes and a 135 bp 59 end fragment of the 23S rRNA (rrlA) gene ( Figure 1). Parts of this region were sequenced from seven additional isolates (accession numbers: HM161644-HM161653; HM194744-HM194759; HM194761-HM194767; HM194760; HM194726-HM194728; HM194731-HM194737; HM194729-HM194730). Alignment of all sequences revealed a low level of polymorphism among isolates in this part of the genome (Figure 1). The arrangement of the genes and their orientation were similar to those found in B. duttonii [16], B. recurrentis [16], B. turicatae and B. hermsii [16,17], and different from the B. burgdorferi species [8,16] in which a 3052-bp region separates the 16S rRNA from the 23S rRNA genes [9,17] (Figure1). We used the sequences of this large contig to generate a phylogenic tree with other TBRF Borrelia for which the equivalent sequences were available. As shown in Figure 2, the Israeli Borrelia isolate is well separated from both Old World and New World Borrelia species and could be defined as a separate species from other TBRF Borrelia species.

Sequencing and analysis of the rrs-ileT intergenic region
The rrs-ileT region sometimes described as intergenic spacer (IGS) has been used in several studies for the taxonomic analysis  Table 1. The phylogenic tree was inferred using the UPGMA method as described in Figure 2. All positions containing gaps and missing data were eliminated from the dataset (Complete deletion option). There were a total of 349 positions in the final dataset. doi:10.1371/journal.pone.0014105.g003 Table 5. Variability and genovar distribution based on the 16S rRNA gene (rrs) sequences of B. persica in Israel and Iran.

Location Genovar
Number of strains 1 Nucleotides at position 2 371 625 1435 1443 The accession numbers of all rrs sequences of Israeli isolates, complete and partial, are listed in materials and methods. Accession numbers of the Iranian strains are U42297 (complete sequence) and EU914141 (partial sequence). 1 The number in parenthesis indicates the number of isolates in each genovar for which the complete rrs sequences are available. 2 Numbering of the nucleotide positions is according to the complete rrs sequence of the R1 strain H1039 (HM161645 Borrelia species [9,18,19,20]. We amplified and sequenced this region in 16 independent B. persica isolates from human blood and from O. tholozani that were collected in Israel and the West Bank. The size of this region was found to be 439 bp in all isolates tested in this work. This is the smallest IGS of known Old World TBRF Borrelia species [18]. The analysis of these sequences revealed three groups with nucleotide substitutions in only two positions ( Table 3). The 439 bp rrs-ileT region contains the alaT gene (Figure 1), but is otherwise non-coding. Surprisingly, in spite of the non-coding nature of most of this region, the level of similarity among the studied strains was high (99.3-100%) whereas with other Borrelia species the level of similarity was considerably lower, ranging between 59.0-70.9% (Table 4). This observation is consistent with the low level of diversity found in IGS sequences of B. hispanica and B. recurrentis [18,19,21], in contrast to the high diversity found in the New World TBRF Borrelia hermsii [22], in Lyme disease Borrelia [9] and in recently deposited IGS sequences obtained (DQ768099-DQ768103, DQ768105) from Israeli human, feline and canine hosts that demonstrated high diversity even among themselves.
Alignment with other Borrelia sequences showed that the IGS locus was especially discriminatory between the local Borrelia isolates and the Old or New World TBRF Borrelia species (Figure 3). The 439 bp rrs-ileT region sequence is not available for the prototype Iranian strains but, as in the case for the entire rrs-rrlA intergenic spacer, the comparison of this sequence with other Borrelia species clearly showed that all Israeli and Palestinian isolates formed a unique cluster separated from all other TBRF Borrelia species (Figure 3).

Sequencing and analysis of the 16S rRNA gene (rrs)
Although analysis of the 7231 bp fragment clearly classifies the HL2610 isolate as an individual species when compared to other TBRF Borrelia species, this complete sequence was not reported for the Iranian reference B. persica strain [7] and therefore could not be used for the definitive classification of the isolate HL2610 as B. persica.
For speciation of the local isolates described in this work, we compared the complete sequence of the rrs gene from the reference B. persica originally isolated in Iran (U42297) [7] and the partial rrs gene sequence of an independent Iranian B. persica strain IRbp1 (U914141) to rrs sequences of local isolates ( Table 1). The complete 16S rRNA genes of three human TBRF Borrelia isolates (HL2610, H1254, H1039) were sequenced ( Table 5). The resulting 16S rRNA sequences were aligned with that of the prototype B. persica Iranian strain (U42297) [7]. Positions of polymorphism between Israeli and the Iranian strain are summarized in Table 3 (C to T at position 1435, A to G at position 1443, a C insertion at position 371 in all the sequences, and in the case of strain H2610 an additional A to G at position 625). The C insertion at position 371 has been also described for both Old and New World TBRF Borrelia species (B. duttonii: U42288, AF107364; B. hispanica: U42294, DQ057988; B. hermsii: U42292, CP000048; B. turicatae: U42299, CP000049). The coefficients of similarity between the local and Iranian strains were very high (99.8%) suggesting that these isolates constitute a single species, different from other TBRF Borrelia which showed coefficients of similarity ranging between 99.1-98.4% (Table 4). As expected, the coefficient of similarity with B. burgdorferi was much lower (96%). These results classify the  Israeli isolates as B. persica. The complete prototype rrs sequence was used by Ras et al. [7] for the classification of B. persica in relationship to other TBRF Borrelia species. Given the low discriminatory power of this conserved gene, the authors concluded that the classification of B. persica as an independent species is relevant but not fully satisfactory [3,7]. Stackebrandt and Ebers suggested that a cutoff range of 98.7-99.0% in the 16S rRNA gene homology is appropriate for species differentiation within a genus [23]. According to these criteria, the data presented here clearly differentiate B. persica from other TBRF Borrelia species and definitively allow the classification of the agent responsible for TBRF in our region as B. persica.
The rrs sequences of the three Israeli isolates, originating from different locations were identical except for one mismatch (A to G at position 625) ( Table 5). Sequence analysis of a 498 bp fragment of the rrs gene comprising the position 625 (amplified by primers 16S1F and 16S7R) from seven additional samples confirmed the existence of two distinct B. persica genovars based on the rrs sequence: genovar R1 with a G and genovar R2 with an A at position 625 (Table 5). Partial sequences of the rrs gene from previously published Israeli isolates (DQ207600-DQ207603, DQ768104, AY763792) showed a similar A to G distribution.
The B. persica strains isolated in Iran, the prototype strain (U42297) and the partial IRbp1 strain (EU914141) have a G at position 625 thus belonging to genovar R1.
The distribution of the local isolates between two genovars based on the A to G modification in the rrs gene places B. persica in a unique phylogenic position. The analysis of the available rrs sequences of TBRF Borrelia indicates that the presence of A or G at this position discriminates between Old World (B. recurrentis, B. duttonii, B. crocidurae and B. hispanica) and New World (B. hermsii, B. parkeri, B. turicatae) TBRF Borrelia, so that a G is characteristic of New World species whereas A is a signature of the Old World TBRF Borrelia. B. persica is the only species were both possibilities are present.
The relationship among the TBRF species is also reflected in the phylogenic tree presented in Figure 4. The analyzed Borrelia species were separated into two clusters, one including the New World species B. turicatae and B. hermsii and, the second dividing into two branches that included on the one hand the African strains B. hispanica, B. duttonii and B. recurrentis and on the other hand all Israeli isolates that grouped together with the reference Iranian B. persica strain.
Sequencing and analysis of the purA gene Among the genes in the large 7231 bp contig, the hypoxanthineguanine phosphoribosyltransferase (hpt), adenylosuccinate lyase (purB) and adenylosuccinate synthase (purA) belong to a group of six open reading frames found in the TBRF Borrelia species but not in B. burgdorferi [24]. The 1284 bp sequence of the purA gene was amplified and sequenced. Analysis of sequences from eight B. persica isolates, seven from human hosts and one from O. tholozani revealed the presence of 4 genovars reflecting five base pair substitutions, two of which result in amino acid substitutions Ser to Leu in genovar P3 and Ser to Pro in genovar P4) ( Table 6). The data obtained from the analysis of purA sequences confirmed that all the local strains belonged to a single phylogenetic cluster ( Figure 5).
The coefficient of similarity between the purA gene sequences of the eight isolates from ticks or human varied from 99.7% to 100% whereas the coefficient of similarity with the purA of other TBRF Borrelia species varied from 83.9% to 89.1% (Table 4) confirming once more the status of B. persica as a distinct TBRF Borrelia species.  Sequence analysis of the flagellin gene (flaB) The sequence of flaB and the sequences of the rrs-rrlA intergenic region are often used for molecular epidemiology of Borrelia ssp. [9,18,25] but unfortunately are not available for the strains from Iran. Analysis of Israeli strains based on the flaB gene has been described before and three distinct genovars (I to III) were defined [2]. Typing is performed by amplifying and sequencing a 750 bp fragment of the flaB gene delineated by primers BOR1 and BOR2 ( Table 2). The sequencing of the partial flaB genes of 12 human and ticks samples originating from Israel and the Palestinian Authority showed that the isolates belonged either to type I or to type II ( Table 1) described previously [2]. The level of similarity (99.0-100%) and modifications found in this study were identical to those described by Assous et al. [2].
Phylogenic analysis based on the partial flaB sequences (data not shown) confirmed that all local isolates from human as well as tick samples formed a unique cluster, different from the New World and Old World clusters [3].

Sequence analysis of the glycerophosphodiester phosphodiesterase (glpQ) gene
Similarly to the purA gene, the presence of the glpQ gene differentiates TBRF Borrelia from Lyme disease Borrelia species that lack these genes [24,26,27,28].
We amplified and sequenced an 813 bp fragment of the glpQ gene from 18 samples isolated in various locations in Israel and the West Bank (Table 1). Polymorphism of the glpQ gene among the 18 samples amplified by PCR is shown in Table 7 with substitutions at position 181 (G to A), 472 (C to T) and 613 (G to T resulting in a change of amino acid Ala to Ser) creating four different genovars (G1 to G4). Recently the sequence of a 668 bp fragment of the glpQ gene from an Iranian strain (IRbp1) was published (accession number EU914143). Although this strain is different from the Iranian reference strain used for sequencing of the 16S rRNA gene, it allows the analysis of phylogenic relationship between Israeli and Iranian Borrelia based on this locus. The coefficient of similarity among the sequences of the Figure 6. Phylogenic tree based on glpQ nucleotide sequences. glpQ sequences of 18 independent isolates from Israel and the West Bank belonging to genovars G1 to G4 were compared to glpQ sequences from other Borrelia species (accession numbers are given in parentheses). The isolates in each genovar are listed in Table 1. The phylogenic tree was inferred using the UPGMA method as described in Figure 2. All positions containing gaps and missing data were eliminated from the dataset (Complete deletion option). There were a total of 637 positions in the final dataset. doi:10.1371/journal.pone.0014105.g006 Figure 7. Phylogenetic tree based on the concatenated alignments of purA, glpQ, flaB and rrs-ileT IGS nucleotide sequences. The tree was inferred using the PHYML program. The percentage of trees in which the associated taxa clustered together in the bootstrap test (100 replicates) is shown next to the branches if higher than 80%. The tree was arbitrarily rooted at midpoint. doi:10.1371/journal.pone.0014105.g007 different isolates from ticks or human whether in the West Bank or in Israel varied from 99.5 to 100%.The of similarity with glpQ sequence from Iran was only 93.8% (Table 4). However, the coefficient of similarity with the glpQ sequences of other TBRF Borrelia species was considerably lower and varied from 89.7% to 86.5% (Table 4). The phylogenic tree generated using these sequences ( Figure 6) showed that isolates in Israel/Palestinian Authority formed a cluster well separated from Old World and New World Borrelia spp. strengthening the classification of the isolates in Israel and in the Palestinian Authority as B. persica. The Iranian strain represented a separate branch, probably reflecting an independent clonal evolution of a subspecies of B. persica; nevertheless this branch was very close to the branch containing all the strains described in this work and the two branches were clearly separated from other Borrelia species (Figure 6).
A partial sequence of the glpQ gene from B. persica isolated in Israel was published before [4]. Its coefficient of similarity to the glpQ sequences in our work is only 91.1%. Surprisingly that glpQ sequence showed a very high degree of similarity (98.4%) with the glpQ gene of B. duttonii with near identity (99.0%) at the protein level. The partial sequence of the rrs gene of this isolate suggests that it belongs to genovar R1. Such a glpQ variant is not represented among the isolates in our collection and may indicate that the diversity among B. persica strains in the Middle East is larger than suggested before.

Multi locus sequence analysis (MLSA)
To definitively confirm the speciation of the Israeli and Palestinian authority Borrelia isolates, MLSA was performed on seven isolates originating from human and tick samples and belonging to both 16S rRNA genovars (Figure 7). The MLSA study was based on genes of different loci both coding (purA, flaB and glpQ) and the non-coding (rrs-ileT IGS) sequences. The alignment of the sequences for each of the loci reveals a low level of diversity among the isolates (Table 4). A concatenated sequence based on these genes was used for a phylogenic analysis with homologous loci of other TBRF Borrelia species. The seven independent local isolates formed a unique cluster separated from all other species (Figure 7).
Sequencing of six or more genes from different loci allows a multilocus sequence analysis (MLSA) that can define bacterial species without sequencing a whole genome or performing DNA-DNA hybridizations [29]. Average nucleotide identity of 95% or more defines a species when using MLSA [29]. Although we used only four loci (besides the 16S rRNA gene), our results show that the average nucleotide identity among the local B. persica isolates was higher than 99% whereas with other Borrelia species, the identity was less than 90% (Table 4), confirming the status of species for B. persica and classifying the Israeli and Palestinian isolates as B. persica independently of the 16S rRNA data presented above.
In this work we were able to classify the human and tick TBRF Borrelia isolates originating from the Palestinian Authority and from Israel as a single defined species B. persica. Indeed, our data confirm for the first time, that B. persica is the infectious agent of TBRF in Israel and the Palestinian Authority. B. persica forms a distinct phylogenetic group that is located between the TBRF Borrelia species of the Old World and the New World. Differences in the flaB, glpQ and IGS loci allowed for more specific typing of the B. persica strains, however there were no noticeable differences between human and tick isolates confirming the co-speciation between B. persica and the vector O. tholozani. Further studies are needed to discover the reservoir of the agent of TBRF in our region.
The data generated in this work are now available for the design of specific primers and probes for comprehensive epidemiological studies and evaluation of diagnostic assays that will improve the detection of the disease in our region.