A glimpse of the bacteriome of Hyalomma dromedarii ticks infesting camels reveals human Helicobacter pylori pathogen.

INTRODUCTION
The tick Hyalomma dromedarii is predominant in camels of Saudi Arabia and harbor multiple pathogens causing disease in humans and animals. Knowing the bacterial community of ticks is crucial for surveillance of known and newly emerging pathogens. Yet, the bacteriome of H. dromedarii remain unexplored to date.


METHODOLOGY
In a cross-sectional survey, we used V3-V4 region of 16S rRNA to characterize the bacteriome of 62 whole H. dromedarii tick samples collected from camels found in Hofuf city in Saudi Arabia.


RESULTS
Sequencing results yielded 217 species incorporated into 114 genera, which in turn belong to the dominant phylum Proteobacteria (98%) followed by Firmicutes (1.38%), Actinobacteria (0.36%), Bacteroidetes (0.17%), meanwhile the phyla Cyanobacteria, Verrucomicrobia and unclassified bacteria were rarely detected. Francisella endosymbiont dominated the bacteriome of H. dromedarii ticks with average abundance of 94.37% and together with Salincoccus sp. accounted for 94.51% of the average sequences. The remaining bacteriome consisted of low abundance of potential pathogens and environmental bacteria. Of these pathogens, we found Helicobacter pylori in the tick H. dromedarii for the first time. Notably, Anaplasma, Ehrlichia and Rickettsia pathogens known to be found in H. dromedarii ticks were not detected.


CONCLUSION
This first preliminary study advances our knowledge about the bacterial community of H. dromedarii ticks and provides a basis for pathogen surveillance and studying the influences of symbionts on vector competence. Presence of pathogens in ticks, raise concerns about potential transmission of these agents to humans or animals.


Introduction
Hyalomma dromedarii is a species of hard-bodied ticks in the family Ixodidae that parasitizes several domestic ungulate animals [1] and the most commonly reported tick attached to camels in Saudi Arabia [2]. The tick H. dromedarii harbors several human and animal pathogens such as Crimean-Congo hemorrhagic fever virus [3], Alkhurma hemorrhagic fever virus [4], Theileria camelensis [5] and Rickettsia species [6]. Thus, H. dromedarii ticks are suspected to play a role in the epidemiology of these pathogens. Notably, most prior studies of H. dromedarii in Saudi Arabia have focused on viral and protozoan pathogens but not on bacterial agents [4][5]. Moreover, globally studies have screened H. dromedarii-borne bacteria but using species-specific PCR-based assay. Consequently, a complete screen of the H. dromedarii bacteriome is needed. Currently, the 16S rRNA metataxonomics analyses circumvent the limitation of previous methods, facilitating detection of more bacterial communities in ticks. From the epidemiological point of view, comprehensive analysis of bacteria residing in ticks of veterinary and medical importance is decisive for monitoring and surveillance of diseases and newly emerging zoonotic pathogens circulating in ticks. So far, the tick microbiome of the genus Hyalomma has only been characterized in the species H. rufipes, H. annotilucm, H. isaaci, H. scupense, H. aegyptium, H. marginatum and H. excavatum [7][8][9] by 16S rRNA metataxonomic approach. Since previous studies already determined that the tick microbiome can alter with a number of factors such as ticks feeding, environment and life stage [10,11], the analysis of these factors is beyond the scope of this study. Therefore, the aim of this study was to survey pathogens of H. dromedarii ticks in the eastern province of Saudi Arabia using 16S rRNA metataxonomic approach.

Tick collection and DNA extraction
In a cross-sectional study, a total of 62 adult H. dromedarii tick isolates were collected from camels at the local animal trade market in Al Hofuf, eastern province, where camels are brought to the market daily from different sites in Saudi Arabia. Samples were collected from April to May in 2017. Ticks were collected from camels and placed in sterile Falcon™ 50mL conical centrifuge tubes. Collected ticks were stored at −20 °C, before whole tick DNA was extracted using the DNeasy Blood &Tissue extraction kit (Qiagen, Hilden, Germany) according to the manufacturer's instructions and kept at −20 °C until used as a template for PCR amplification.

Tick identification
Morphological and molecular identification were used for tick identification. PCR amplifications of ribosomal 16S rRNA were generated with published primers [12]. Two microliters of tick DNA 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, Subang Jaya, Malaysia), 5 μL of 10X ViBuffer (Vivantis Technologies, Subang Jaya, Malaysia) and 2 μL of dNTPs (10 mM). The volume was adjusted to 25 μL by adding distilled water. Thermal cycling was performed on a Tpersonal Thermocycler (BIOMETRA, Gottingen, Germany) with an initial 15-minutes cycle at 95°C followed by 35 cycles consisting of 30 seconds at 94°C, 1 minute at 55°C or 60°C depending on the primer and 1 minute at 72°C, followed by a 10 minutes final extension step at 72°C. To rule out DNA or amplicon contamination, molecular grade water negative control was used throughout the steps of the protocol. . Ticks were grouped based on sex and engorged or non-engorged ticks. Each pool composed of male and female plus engorged and nonengorged ticks. The number of males, female, engorged and non-engorged ticks between pools was not the same due to their unequal number in the collected samples. For instance, the 62 adult ticks consisted of 27 engorged and 35 non-engorged ticks. Therefore, we have unequally number of engorged and non-engorged ticks. This distribution will not affect our main goal which is to screen pathogen and not to estimate difference between female and male ticks or engorged and nonengorged ticks which is published previously in other tick species.
Briefly, The V3-V4 segment of the 16S rRNA gene was amplified with Bakt_341F and Bakt_805R primers [14]. The amplicon library was constructed by ligating sequencing adapters and indices to purified PCR products using the Nextera XT Kit (Illumina, San Diego, CA, USA) according to the 16S rRNA metataxonomics sequencing library preparation protocol (Illumina, San Diego, CA, USA). Then libraries concentration was measured, equimolar volume of each of the libraries was pooled and sequenced on an Illumina's Miseq platform with pairedend 300 bp reads by Macrogen Inc (Seoul, South Korea). MiSeq reads were assembled by FLASH version 1.2.11 [15] which merge overlapping pairedend reads. Read trimming, filtering with a quality score offset 33 and OTU picking with a 97% identity cut-off was performed using CD-HIT-OTU software [16]. OTUs were classified by blast against NCBI 16S rRNA database with BLASTN using default parameters [17]. QIIME software was used to assign taxonomy and perform rarefaction curves and alpha diversity analyses (Chao1 index and sample coverage) [18].
For species-level identification using V3-V4 16S rRNA sequences region, Villmones et al. 2018 [19] recommends ≥ 99.3% similarity with a trusted reference species together with a minimum distance of > 0.8% to the closest species. Based on the levels of intra-species sequence variation we observed in Genbank sequences, we adopted a more stringent cut off ≥ 1 % minimum distance to the closest species while keeping a similarity of ≥ 99.3%.

Results
Taxonomic analysis of 16S rRNA sequencing data All tick samples collected in this study were genetically identified by sequencing of partial 16S rRNA gene as H. dromedarii ticks. For 62 H. dromedarii tick DNA samples (2 individual, 8 pooled), we obtained after removal of low quality and chimeric reads, a total of 755,940 high quality reads. The observed rarefaction curves and chao1 rarefaction curves reached plateau for all samples (Figure 1). The Good's coverage estimates range between 0.99% to 1.00%. These results show that the sequencing depth was sufficient to estimate 99% of the bacterial diversity and species richness in all samples ( Figure 1).
A total of 546 OTUs were identified at 97% sequence similarity level, which were assigned to 6 phyla, 70 families and 114 genera. The unclassified OTUs at the phylum and the genus level were 12 and 29 OTUs respectively ( Table 1). The bacterial calculated richness varied from (51 to 57 OTUs) per individual samples and (26 to 106 OTUs) per pooled samples ( Table 1). The lowest number of observed OTUs was 26 in sample 55, whereas the highest was 106 in sample 30 (Table 1).
At the phylum level, Proteobacteria was found to be the most dominant with average abundance of 98.12% followed by Firmicutes (1.34%), Actinobacteria (0.33%), Bacteroidetes (0.16%), meanwhile the phyla cyanobacteria, Verrucomicrobia and unclassified were rarely detected ( Table 2). Among these phyla, Proteobacteria, Firmicutes and Actinobacteria were found to be present in all samples. The total number of bacteria assigned to Firmicutes was 214 OTUs followed by Proteobacteria (160 OTUs) and Actinobacteria (99 OTUs) other phyla contain less than 50 OTUs were shown in (Table 1).   Table 1). Other genera having abundance less than 0,1% were listed in (Supplementary Tables 2 and 3).
At the species level, thirty-three out of 114 genera contained more than one species, of these, Corynebacterium, Bifidobacterium and Bacillus genera constituted the most diverse genera, each containing 11 species. Other genera were listed in (Supplementary  table 1, 2 and 3). Following the classification criteria adopted in the study, only 16 out of 217 unique OTUs could be classified to the species level (Table 3). Of these, 2 out of 3 OTUs of the genus Helicobacter were classified as Helicobacter pylori with similarity value of 99.7% to H. pylori in GenBank database and had a distance of more than 1 to Helicobacter cetorum and H. pullorum with 98.1% identity. The V3-V4 16S rRNA sequence phylogeny clustered H. pylori sequence found in the study with Genbank H. pylori sequences, while separated from other Helicobacter species ( Unfortunately, the Francisella V3-V4 region of 16S rRNA gene matches best to Francisella-like endosymbiont (FLE) of Hyalomma marginatum at 98.71% identity, which is below the threshold 99.3% for species assignment adopted in this study. To get better taxonomic resolution a 1071bp region of and Francisella hispaniensis (CP018093.1) respectively. As for Francisella tularensis strains, the similarity varied from 98.11% to 97.9%. The phylogenetic inference based on 1071bp 16S rRNA gene sequence showed that Francisella sp. was genetically related to others symbiotic Francisella, while separated from the others pathogenic Francisella species (Figure 3, Supplementary Figure 2).   Table  2). It accounted for 1.15% of average sequences and 24.42% (53 species) of total bacteria and 3) low prevalent bacterial genera having prevalence less than 40%, which represent 3.5% of average sequences and constituted 66.36% (144 species) of total bacterial (Supplementary Table 3).

Discussion
In this study, we conducted a cross sectional survey for bacterial community in whole H. dromedarii ticks in Saudi Arabia. The bacteriome of H. dromedarii ticks was analyzed via sequencing of the V3-V4 segment of 16S rRNA gene using Illumina MiSeq sequencer. Although ticks were randomly collected from camel, H. dromedarii ticks are the only ticks found during collection. This finding is consistent with previous report revealing the dominance of H. dromedarii ticks in camels of Saudi Arabia [3].
Using this approach, a total of 6 phyla were detected in tick H. dromedarii and Proteobacteria was the most abundant which agreed with the composition of bacterial community reported in several tick species [20]. Although Proteobacteria is more dominant at sequences level than Firmicutes phylum, the number of genera assigned to Firmicutes (214 OTUs) surpassed the number of genera assigned to Proteobacteria (160 OTUs). Furthermore, Proteobacteria, Firmicutes and Actinobacteria phyla were prevalent in all samples either in individual or pooled samples.
Several human pathogens were detected in this study such as Helicobacter pylori, the causative agents of stomach peptic ulcer disease in human. A disease proposed to spread among human through the oral-oral or fecal-oral routes [21]. H. pylori were detected in 3 tick pools (30%). However, the presence of H. pylori in H. dromedarii ticks does not prove that the ticks act as reservoir or competent vector, but remains to be elucidated. It also raises a question how H. dromedarii ticks acquire this bacterium. Previous studies reported H. pylori in the stomach of domestic animals without having gastritis [22,23] and also found in the milk of camel, cow and goat. Therefore, proper caution is required when removing or handling ticks during collection to avoids hand contamination. To our knowledge, this is first report of H. pylori DNA in ticks, but non-H. pylori species such as Helicobacter bizzozeronii was detected before in H. rufipes ticks in china [7].   (Table 4). Notably, Staphylococcus and Corynebacterium have previously been detected from the saliva content of H. flava [24]. Among moderately prevalent genera we also found Pseudomonas and Enterococcus which were detected previously from internal organ of Rhipicephalus (Boophilus) microplus, Ixodes ovatus and I. persulcatus [25]. Thus, these genera are probably true members of tick bacteriome. Furthermore, the species Acinetobacter variabilis, Acinetobacter schindleri, Corynebacterium confusum and Corynebacterium massiliense have been detected before in human clinical specimens [26,27]. Among low prevalent genera we also found pathogenic genera such as Klebsiella, Lactococcus, Lysinibacillus and Massilia that have species previously been detected in human clinical specimens [28][29][30]. The lower prevalence and abundance of pathogenic genera suggest they are likely transient bacteria acquired from surrounding environment. Other bacterial species were commensal of human intestinal tract, such as, Granulicatella adiacens, Akkermansia muciniphila and anaerostipes hadrus that is found in human colon (Table 3

) [31-33].
Coxiella burnetii is the etiological agent for Q fever, a worldwide zoonotic disease reported in human and animals such as camels, sheep, goats and cattle [34]. Transmission of Q fever to animals via tick bite in nature has not been confirmed, yet ticks have been experimental shown to be competent vectors for the transmission of C. burnetii to animal hosts [35]. However, the presence of C. burnetii in H. dromedarii and some other ticks suggest a role for ticks in the epidemiology Q fever. Most noticeable in our survey is the absence of the genus Coxiella in our bacteriome analysis although the infection is common in camels in Saudi Arabia [36]. The other noteworthy tick-borne pathogens found previously in H. dromedarii ticks outside Saudi Arabia but not found in our bacteriome analysis include Rickettsia, Anaplasma and Ehrlichia genera [6,37]. It is possible that these pathogens are absent in H. dromedarii tick population in Saudi Arabia, or have low prevalence, thus not detected here because of the small sample size. In addition to pathogenic and endosymbiotic bacteria, environment associated bacteria were observed including the soil bacteria Solibacillus which was isolated previously from the midgut of sand flies [40] and salt mine associated Salinicoccus kunmingensis. Other soil members among low prevalent genera comprise the genera Pusillimonas, Oxalicibacterium and N. islandensis [41]. Furthermore, S. mesophilum and D albogilva and L. defluvii [42][43][44] were reported from wastewater. Although the detected genera Lysobacter, Pantoea, Paracoccus, Pontibacter and Pseudomonas are frequent members of sandy soil of Saudi Arabia, the current study has yet to determine that these genera are acquired from the environment [45]. Finally, the presence of 12 unclassified OTUs may indicate the existence of as yet uncharacterized novel species.
The high prevalent bacterial genera (Francisella, Salinicoccus, Corynebacterium and Staphylococcus) probably encodes certain functions associated with tick survival and reproduction, which warrants further investigation to elucidate. Tick as external parasite can acquire bacteria from host skin or environment as shown previously for H. dromedarii ticks [46]. Hence, environmental factors may shape H. dromedarii bacteriome. In sum, our review of 12-tick species microbiota reveals that the genera Bacillus, Staphylococcus Corynebacterium and Pseudomonas are highly prevalent in ticks. Furthermore, 46.5 % (53 genera) of bacteria found in the present study have been detected previously from other tick species (Table 4).
As for bacterial abundance, FLE of H. dromedarii was the most abundant coexisting bacteria present in all samples, with abundance ranged from 93% to 99% of sequence reads in pooled samples except sample 55 having abundance of 70%, while in individual samples the abundance was 95.6 and 96.6% in, Other genera with abundance above 1% include Proteus (29.7%), Acinetobacter (3.38%) and Staphylococcus (2.22%), while the remainder of detected bacterial genera had abundance less than 1%. The noteworthy in sample 55 it exhibits the lowest species richness (26), high abundance of Proteus (29.70%) and low abundance of FLE of H. dromedarii (70.03%), contrary to high abundance of FLE of H. dromedarii (> 93%) in the rest of samples. Although this observation warrants some sort of correlation analysis of absolute abundances across individuals to explain this finding, previous studies have demonstrated that tick bacteriome can interfere with pathogens colonization and transmission. For instance, in Dermacentor andersoni ticks, a reduction in Francisella endosymbionts was associated with lower Francisella novicida abundance levels [47]. In addition, an increase in ratio of Rickettsia bellii was associated with reduction of Anaplasma marginale levels in Dermacentor andersoni ticks [47]. Furthermore, a study reported that Ixodes scapularis microbiome composition could influence Borrelia burgdorferi colonization [48].
One of the limitations of our study it focused on whole tick bacteriome, therefore it did not differentiate between internal bacteria of ticks and the bacterial species residing on the exoskeleton. Regardless of this, 35.1% of genera detected in the current study have been reported from internal organs (saliva, midgut and ovaries) of several ticks (Table 4). On the other hand, biologically transmission is not the sole route of pathogen transmission; non-salivary mechanical bacterial transmission can also occur by contamination of injuries induced at feeding site with exoskeleton bacteria, raising the importance of exoskeleton bacteria. Another limitation of our study is small sample size, which prevents us from confirming the presence or absence of some tick-associated pathogens such as Rickettsia, Anaplasma and Ehrlichia.

Conclusion
This study has characterized the bacteriome of whole H. dromedarii ticks and revealed that the ticks mainly colonized by Francisella endosymbionts with low abundance of potential pathogens and environmental bacteria. However, pathogens detected herein do not indicate that H. dromedarii is a competent vector but may pose potential risk to humans and animals.

Annex -Supplementary Items
Supplementary Figure 1. Neighbor-Joining tree with branch lengths of the V3-V4 16S rRNA sequences of Helicobacter species and OTU identified as Helicobacter pylori (Green circle) in this study. The bootstrap consensus tree inferred from 100 replicates. Bootstrap values > 50% are shown. The tree is drawn to scale, with branch lengths in the same units as those of the evolutionary distances used to infer the phylogenetic tree. The evolutionary distances were computed using the Tamura-Nei method and are in the units of the number of base substitutions per site. The analysis involved 44 nucleotide sequences. All positions containing gaps and missing data were eliminated. There were a total of 431 positions in the final dataset. The analyses were conducted in MEGA7 [58].

Supplementary Figure 2.
Maximum-likelihood tree with branch lengths based on 1071bp 16S rRNA gene sequences of FLE Hyalomma dromedarii (sample no 52 and 56) generated as part of this study and selected sequences of Francisella species from the GenBank. The tree is drawn to scale, with branch lengths measured in the number of substitutions per site (next to the branches). The analysis involved 36 nucleotide sequences. There were a total of 1070 positions in the final dataset. Analysis was conducted in MEGA7 [58] based on Tamura-Nei model with 100 bootstraps.