Low genetic diversity among Francisella-like endosymbionts within different genotypes of Hyalomma dromedarii ticks infesting camels in Saudi Arabia

Background and Aim: Hyalomma dromedarii ticks are vectors of disease agents and hosts of Francisella-like endosymbionts (FLEs). Knowledge about intraspecific genetic variation among H. dromedarii and its Francisella species is limited. The aims of this study were to investigate whether certain H. dromedarii genotypes are specialized in carrying specific Francisella species genotypes and scrutinize the population structure of H. dromedarii ticks in Saudi Arabia. Materials and Methods: We collected 151 H. dromedarii ticks from 33 camels from 13 locations in Saudi Arabia. The second internal transcribed spacer (ITS2), cytochrome c oxidase subunit-1(COI), and 16S rRNA genes were used for single- and multi-locus sequence typing and phylogenetic analyses. H. dromedarii-borne Francisella was screened using the tul4 gene and 16S rRNA Francisella-specific primers followed by amplicon Sanger sequencing. Results: Single-locus typing of ticks using ITS2, 16S rRNA, and COI genes yielded 1, 10, and 31 sequence types (ST), respectively, with pairwise sequence similarity of 100% for ITS2, 99.18-99.86% for COI, and 99.50-99.75% for 16S rRNA. COI sequence analysis indicated a lack of strict geographical structuration, as ST15 was found in both Saudi Arabia and Kenya. In contrast, multilocus sequence typing resolved 148 H. dromedarii ticks into 39 genotypes of ticks and three genotypes of FLEs. The ST2-FLE genotype was carried by the tick genotype ST35, while the ST1-FLE genotype and 41.89% of the ST3-FLE genotype were carried by the tick genotype ST32. Accordingly, there appeared to be no specialization of certain tick genotypes to harbor-specific FLE genotypes. Conclusion: For the 1st time, we have provided an overview of the population structure of H. dromedarii ticks and FLE strains. We found a low level of genetic diversity among FLEs and non-specialized circulation of FLEs among H. dromedarii ticks.


Introduction
Hyalomma dromedarii are ticks that infest camels and a wide range of other ungulate animals, including cattle and sheep [1,2]. The geographical distribution of H. dromedarii extends from North, Northwest, Central, and East Africa to the Middle East and Central and South Asia [3]. It is the most commonly reported tick seen attached to camels in Saudi Arabia [4]. H. dromedarii harbors several bacterial, viral, and parasitic pathogens that can infect humans and other animals. For example, in Saudi Arabia, H. dromedarii harbors Sindbis virus, Chick Ross virus, Kadam virus, and Alkhurma hemorrhagic fever virus; it is, therefore, suspected to play a role in the epidemiology of these pathogens [5]. In Iran, [6] molecular evidence suggests that Crimean-Congo hemorrhagic fever virus is present in H. dromedarii. Recently, in Pakistan [7], Theileria annulata parasites were found in H. dromedarii ticks. Notably, it is anticipated that new pathogens will emerge. As well as pathogenic microbes, H. dromedarii carries non-pathogenic microbes and probably endosymbiont organisms, such as Francisella-like endosymbionts (FLEs) [8,9]. A recent study of H. aegyptium ticks showed that Francisella spp. were significantly more abundant than other bacteria [10]. Furthermore, many symbiotic bacteria in ticks are known to interact with other microbes, affecting the competence of vectors and the ability of microbes to colonize and persist within ticks [11,12]. Hence, investigating the prevalence and diversity of FLEs in ticks is of great interest.
Intraspecies genetic analyses of H. dromedarii are important for defining its population structure. Different tick genotypes may differ in their pathogen Available at www.veterinaryworld.org/Vol.13/July-2020/29.pdf colonization and competence. The veterinary and medical importance of this vector emphasizes the need for global phylogenetic studies to advance our knowledge of the relationship between different tick genotypes and their microbiota. To date, at the global level, few molecular studies have been performed to examine the intraspecies genetic diversity of H. dromedarii; those studies that have been performed were limited either by small sample sizes or the number of loci tested. One of these studies, in Kenya, revealed the presence of cryptic hybridization in H. dromedarii ticks [13]. In another study, H. dromedarii samples from Egypt were genotyped by sequence analysis of 18S rDNA, the cytochrome c oxidase subunit-1 (COI) gene, and 16S rDNA [14,15]. In India, H. dromedarii ticks were identified by the COI gene; further characterization and establishment of their phylogenetic status were achieved by sequencing the calreticulin gene and the second internal transcribed spacer (ITS2) [16]. Despite the importance of Francisella species and H. dromedarii ticks, there is limited information about intraspecies genetic diversity among H. dromedarii and their associated Francisella species globally and in Saudi Arabia specifically.
In this study, we attempted to address these issues by performing sequence analysis of the Francisella 16S rDNA gene, which has been shown to be informative for population analyses of FLEs [17]. For ticks, we selected the most commonly used genetic markers (the COI gene, ITS2, and 16S rDNA sequences) to facilitate global comparisons between geographically distant regions. The aims of this study were to investigate whether certain H. dromedarii genotypes are specialized in carrying specific Francisella species genotypes and scrutinize the population structure of H. dromedarii ticks in Saudi Arabia.

Ethical approval
Ethical approval is not necessary for this study.

Tick collection and DNA preparation
One hundred and fifty-one adult ticks, obtained from 33 camels (20 males and 13 females), were collected from 13 locations in Saudi Arabia between April 2017 and January 2020. The locations are indicated on the map shown in Figure-1. It was not possible to collect the same number of samples from each location as sampling was dependent on the availability of ticks in the region. Ticks were removed from camels by hand and stored in sterile 50 ml Falcon™ conical centrifuge tubes. Ticks were placed in separate tubes according to their attachment site on a camel. The tubes were stored at −20°C until whole tick DNA was extracted using a DNeasy Blood and Tissue Kit (Qiagen, Hilden, Germany), according to the manufacturer's instructions. Extracted DNA was stored at −20°C until being used for polymerase chain reaction (PCR) amplification and sequencing.

Tick genotyping
PCR amplification of ITS2 of the rRNA gene and two mitochondrial genes, COI and 16S rRNA, was performed, using both specifically designed primers and previously published primers (Table-1) [18][19][20]. H. dromedarii DNA (2 μl) and 0.3 μl of each primer (10 pmol) (Macrogen, Seoul, South Korea) were added to the PCR mixture, containing one unit of Max Taq DNA Polymerase (Vivantis Technologies, Malaysia), 5 μl 10× ViBuffer (Vivantis Technologies, Malaysia), and 2 μl dNTPs (10 mM). The volume was adjusted to 25 μl by adding distilled water. Thermal cycling was performed on a TPersonal Thermocycler (BIOMETRA, Germany) with an initial 15 min cycle at 95°C, followed by 35 cycles consisting of 30 s at 94°C, 1 min at 55°C or 60°C depending on the primers, and 1 min at 72°C, followed by a final 10 min extension step at 72°C. To rule out DNA or amplicon contamination, Molecular Grade Water was used as a negative control throughout each step of the protocol. The PCR amplicons were sequenced by Macrogen Inc. (Seoul, South Korea) using BigDye (Applied Biosystems, California, USA) on an ABI3730XL  Available at www.veterinaryworld.org/Vol.13/July-2020/29.pdf DNA sequencer (Applied Biosystems, California, USA).

Francisella detection and genotyping
The detection of Francisella DNA in individual ticks was performed through PCR of the total DNA of individual ticks, using the 16S rRNA Francisella-specific primers Fr153F0.1 (5-GCC CATTTGAGGGGGATACC-3) and Fr1281R0.1 (5-GGACTAAGAGTACCTTTTTGAGT-3) [21], as follows: H. dromedarii DNA (3 μl) and 0.3 μl of each primer (10 pmol) (Macrogen Inc., South Korea) were added to the PCR mixture, which contained one unit of Max Taq DNA Polymerase (Vivantis Technologies, Malaysia), 5 μl 10X ViBuffer (Vivantis Technologies, Malaysia), and 2 μl dNTPs (10 mM). The volume was adjusted to 25 μl by adding distilled water. Thermal cycling was performed on a TPersonal Thermocycler (BIOMETRA, Germany) with an initial 15 min cycle at 95°C, followed by 35 cycles consisting of 30 s at 94°C, 1 min at 55°C, and 1 min at 72°C, followed by a final 10 min extension step at 72°C. Molecular Grade Water was used as a negative control throughout each step of the protocol to rule out DNA or amplicon contamination. PCR product purification and sequencing were performed by Macrogen Inc. using BigDye (Applied Biosystems, California, USA) on an ABI3730XL DNA sequencer (Applied Biosystems, California, USA). The obtained sequences were blasted against the NCBI non-redundant (nr) database to identify the species they were closest to.
PCR amplification of Francisella 17 kDa membrane lipoprotein (tul4) was generated using the primers designed herein: tul4HdF (5-CTAATCCCGAAATAATATTGATAGGT-3) and tul4HdR (5-CAGTTGCCCAAGTCTTATCATTC-3). The PCR parameters used were the same as those used for the 16S rRNA Francisella PCR protocol used in the current study.

Data collection and sequence analysis
Sequences of tul4, ITS2, COI, and 16S rRNA of H. dromedarii were retrieved from GenBank and compared with the sequences obtained in the present study. The nucleotide sequences were assembled and edited in CLC Main Workbench 8 software (CLC bio, Aarhus, Denmark). Sequences of the same gene were aligned using the algorithm in CLC Main Workbench 8 with default parameters. All sequences of the same gene were trimmed to the same length to enable comparisons. The ITS2, COI, and 16S rRNA nucleotide sequences analyzed herein were concatenated. Each particular sequence of COI, ITS2, and the 16S rRNA gene, in addition to the concatenated sequence, was assigned a sequence type (ST) number. The marker discrimination power of the DNA markers, COI, ITS2, and the16S rRNA gene, was calculated using the Hunter-Gaston index [22].
Phylogenetic reconstruction was performed separately for the ITS2, COI, and 16S rRNA gene sequences. In addition, relationships among H. dromedarii ticks sequenced herein and H. dromedarii ITS2, COI, and 16S rRNA sequences of H. dromedarii ticks obtained from GenBank were inferred using Bayesian phylogenetic analyses in MrBayes v3.2.6, as follows: Nucleotide sequences were first aligned using MUSCLE software. Bayesian phylogenetic analyses were then performed using MrBayes v3.2.6 [23] with a Jukes-Cantor model of nucleotide substitution. Sampling was performed using three independent runs (each having one cold chain and three heated chains), which were run for 1,000,000 generations.

Nucleotide sequence accession numbers
All sequences are available through GenBank. For the COI gene sequence, the accession numbers are MH590861 to MH590886; for 16S RNA, the accession numbers are MH569476 to MH569481; and for ITS2, the accession number is MH571755.

Hyalomma species
Of 151 ticks (engorged and non-engorged) collected from 33 camels (20 males and 13 females), 148 (98.01%) were H. dromedarii, 2 (1.32%) were H. impeltatum, and 1 (0.66%) was Hyalomma spp. Tick occurrences per camel were rare due to the use of acaricides. The ticks were mainly found attached to the brisket, axilla, lower eyelid, jaw, base of tail, and perineum region. The two H. impeltatum ticks were removed from the base of the tail of a female camel in Afif ( Figure-1), while the single Hyalomma spp. was removed from the base of the tail of a female camel in Riyadh (Figure-1). H. dromedarii ticks were recovered from all attachment sites mentioned above.

Genotyping of ticks
For all PCR experiments, the 148 DNA samples were efficiently amplified by PCR using the three genetic markers, and the negative controls remained negative. The lengths of the assembled sequences after trimming low-quality sequences were 722 bp for CO1 and 401 bp for 16Sr RNA, making it easy to obtain accurate sequences using the primer pairs COX-F/ COX-R and 16SrRNA-F/16SrRNA, respectively. However, for ITS2, we failed to assemble high-quality sequences using ITS2-F/ITS-2R primers only. Specific internal primers (480-F, 480-R, and 720-F) were, therefore, designed, enabling assembly of a 1340 bp sequence.   (Table-4). The genotype FLE-H. dromedarii ST-3 represented 98.6% (146 out of 148 samples), followed by genotype FLE-H. dromedarii ST1 (0.7%), which was found in the Buraidah region, and FLE-H. dromedarii ST2 (0.7%), which was found in the Asir region (Figure-1).

Phylogenetic analysis
Consultation of the GenBank database revealed that H. dromedarii sequences comprised 14 16S rRNA nucleotide sequences, from Senegal, Pakistan, Egypt, Saudi Arabia, and Iraq. There were 22 COI nucleotide sequences, from Kenya, India, United Arab Emirates, Ethiopia, Senegal, Pakistan, Egypt, Saudi Arabia, and Iraq. ITS2 comprised 15 nucleotide sequences, from India, Egypt, Pakistan, Iraq, Saudi Arabia, Iran, and Iraq. Not all sequences retrieved from GenBank were included in this analysis; short sequences were excluded. The concatenation of the 16S rDNA and COI sequences of H. dromedarii separated the strains into different clades, as shown in Figure-2.
A phylogenetic tree combining the STs of 16S rRNA sequences obtained in this study with GenBank 16S rRNA sequences yielded three clades: Clades A and B, consisting of isolates from Egypt and clade C, consisting of Saudi Arabia H. dromedarii isolates clustered with ticks from Senegal, Iraq, and Pakistan (Figure-3). The phylogenetic analysis combining STs of COI sequences found in ticks in our study with GenBank COI sequences revealed an absence of strict geographical structuration since genotype ST15 was found in both Saudi Arabia and Kenya (Figure-4).
For FLEs of H. dromedarii, the phylogenetic tree inferred using the tul4 sequences clearly separated FLEs of Hyalomma strains from the FLEs of other ticks (Figure-5).

Discussion
The multityping approach used in this study revealed a high degree of genetic diversity among H. dromedarii collected from camels in Saudi Arabia. Even though we examined a relatively small set of H. dromedarii, we observed a high diversity index among these ticks, since 148 H. dromedarii ticks yielded 39 genotypes, as determined by the concatenation of COI and 16S rRNA gene sequences. The use of several genetic markers offers better resolution over the use of a single genetic marker, as has been used previously [14,15]. Our results revealed that COI sequences were more variable, in contrast with the low nucleotide variation seen in 16S rRNA  and the highly conserved ITS2 sequence. These findings were consistent with those of previous tick studies [24][25][26]. Consequently, we concluded that COI is a more informative marker than the 16S rRNA gene for studying intraspecific diversity, whereas the ITS2 sequence is not informative but is suitable for H. dromedarii identification. The tick genotype ST32 represented the dominant genotype prevalent in 11 out of 13 locations, which suggests the circulation of this genotype throughout Saudi Arabia. Plausible explanations for this could be that the ST32 genotype is well adapted to persist in its camel host or the frequent movement of camels between different locations as they seek pasture. Furthermore, ST32 showed no preference for a specific attachment site and was found widely distributed around a camel's body.
A 16S rRNA phylogenetic tree showed that H. dromedarii ticks of Saudi Arabia were clustered with ticks from Senegal, Iraq, and Pakistan. Moreover, COI phylogenetic analysis revealed an absence of strict geographical structuration, as genotype ST15 was found in both Saudi Arabia and Kenya. A reasonable explanation for this could be the longstanding trade in camels between these countries, which might play a significant role in spreading ticks. More H. dromedarii samples are, therefore, needed to shed light on the global diversity of this species.
In this study, a total of 148 ticks were screened for Francisella species. F. tularensis was not detected in any samples. In contrast, FLEs were highly prevalent (100%) in H. dromedarii. These results indicate some sort of beneficial relationship between FLEs and H. dromedarii. The previous studies into H. dromedarii showed differences in the prevalence rates of FLEs, which varied from 6% to 89.8% [8,9]. A possible explanation for this discrepancy in results could be explained either by FLEs being present but undetectable due to low concentrations of DNA or that certain genotypes of H. dromedarii ticks have no capacity to harbor FLEs.
The presence of different STs of FLEs of D. andersoni and D. variabilis has been reported in Canada following the sequencing of 16S rRNA [17]. In addition, the FLE-ST2 genotype has been reported in Egypt. Here, we observed no intraspecies discrimination among Saudi Arabian FLE strains using 16S rRNA genes, in contrast with the tul4 gene, which grouped the 148 FLEs of H. dromedarii into three genotypes.  Available at www.veterinaryworld.org/Vol.13/July-2020/29.pdf  Available at www.veterinaryworld.org/Vol.13/July-2020/29.pdf