Co-Occurrence of Francisella, Spotted Fever Group Rickettsia, and Midichloria in Avian-Associated Hyalomma rufipes

The migratory behavior of wild birds contributes to the geographical spread of ticks and their microorganisms. In this study, we aimed to investigate the dispersal and co-occurrence of Francisella and spotted fever group Rickettsia (SFGR) in ticks infesting birds migrating northward in the African-Western Palaearctic region (AWPR). Birds were trapped with mist nests across the Mediterranean basin during the 2014 and 2015 spring migration. In total, 575 ticks were collected from 244 birds. We screened the ticks for the species Francisella tularensis, the genus Francisella, and SFGR by microfluidic real-time PCR. Confirmatory analyses and metagenomic sequencing were performed on tick samples that putatively tested positive for F. tularensis during initial screenings. Hyalomma rufipes was the most common tick species and had a high prevalence of Francisella, including co-occurrence of Francisella and SFGR. Metagenomic analysis of total DNA extracted from two H. rufipes confirmed the presence of Francisella, Rickettsia, and Midichloria. Average nucleotide identity and phylogenetic inference indicated the highest identity of the metagenome-assembled genomes to a Francisella-like endosymbiont (FLE), Rickettsia aeschlimannii, and Midichloria mitochondrii. The results of this study suggest that (i) FLE- and SFGR-containing ticks are dispersed by northbound migratory birds in the AWPR, (ii) H. rufipes likely is not involved in transmission of F. tularensis in the AWPR, and (iii) a dual endosymbiosis of FLEs and Midichloria may support some of the nutritional requirements of H. rufipes.


Introduction
Ticks (Acari: Ixodida) transmit pathogens of both human and veterinary importance, such as bacteria in the genera Anaplasma, Borrelia, Coxiella, Francisella, and Rickettsia. They also can be co-infected with different pathogens that potentially can cause co-infections in hosts [1,2]. Additionally, ticks harbor endosymbionts living symbiotically within them. Bacterial endosymbionts of ticks are mostly from the genera Coxiella, Francisella, and Rickettsia; they are closely related to pathogens and may be necessary for the survival of the host [3]. Ticks are strictly hematophagous, meaning their diet consists solely of vertebrate blood, which is nutritionally unbalanced since it contains a high level of proteins but few vitamins [4]. Endosymbiotic bacteria present within tick cells are believed to support the dietary requirements of ticks by providing nutrients that are absent in vertebrate blood [5,6].
The bacterium Midichloria mitochondrii was first described as an endosymbiont of Ixodes ricinus ticks known to reside primarily in the ovarian primordia or ovaries, enter the mitochondria, and be transmitted by transovarial transmission to the offspring [26][27][28]. The prevalence of M. mitochondrii in female I. ricinus ticks has been reported to be 100% [27]. The bacterium also has been detected in ticks of the genera Hyalomma, Rhipicephalus, Amblyomma, and Haemaphysalis [29]. DNA of M. mitochondrii and of bacteria related to M. mitochondrii has been detected in the salivary glands of I. ricinus ticks [30] and in blood samples from canines [31], respectively, and antibodies against Midichloria have been detected in blood samples collected from canines and humans bitten by ticks [31,32], indicating potential horizontal transmission of the bacterium.
The genus Rickettsia has been divided into four groups [33]. Most tick-borne Rickettsia belong to the spotted fever group (SFG), which includes members that are considered to be emerging human pathogens in Europe (e.g., Rickettsia conorii, Rickettsia massiliae, and Rickettsia aeschlimannii) [34,35]. Rickettsia species in the SFG are transmitted to humans by multiple tick genera, including Rhipicephalus, Ixodes, and Hyalomma [34]. Wild birds are frequently parasitized by Ixodes and Hyalomma ticks, and the migratory behavior of the avian hosts aids in the geographical spread of ticks and their associated microorganisms [36][37][38]. Because they are intracellular tick-borne bacteria, Francisella, Midichloria, and Rickettsia are difficult to culture, and culture-independent generation of genome sequences is of importance for increasing the knowledge and understanding of these bacteria. In this study, we aimed to investigate the dispersal and co-occurrence of Francisella and SFG Rickettsia (SFGR) species in ticks infesting northbound migrating birds in the African-Western Palaearctic region (AWPR), using microfluidic real-time (q) PCR and metagenomics.

Trapping of Birds and Collection of Ticks
Birds were trapped using mist nets at bird observatories in Spain (several sites in the provinces of Huelva and Sevilla and the Canary Islands: 37 • [39] for additional details.

DNA Extraction
In brief, absolute ethanol (Sigma-Aldrich, Merck, Darmstadt, Germany) and sterile H 2 O were used for surface sterilization of the ticks before homogenization. Mechanical homogenization was performed using a stainless-steel bead (Qiagen, Hilden, Germany), TRIzol TM (Invitrogen, ThermoFisher Scientific, Waltham, MA, USA), and a TissueLyser II (Qiagen). After homogenization, additional TRIzol TM was added to the homogenate, followed by centrifugation and collection of the supernatant. RNA was isolated and removed using a phase separation technique, in which chloroform (Sigma-Aldrich) was added to the supernatant. Thereafter, DNA was extracted from the organic phase using a back-extraction buffer, inversion, and centrifugation. DNA present in the upper phase was purified using the Nucleospin gDNA Clean up kit (Macherey-Nagel, Bethlehem, PA, USA). The DNA was eluted in DE buffer and stored at −20 • C. For additional details, see Hoffman et al. 2021 [39].

Francisella
Tick extracts were screened for the presence of DNA from the genus Francisella and the species F. tularensis specifically by microfluidic qPCR (BioMark TM Dynamic Arrays, Fluidigm, CA, USA) and with FopA (genus-specific) and Tul4 (species-specific) primers and probes (Table 1), at the Animal Health Laboratory (Paris, France) according to Michelet et al. [40]. Samples with a cycle threshold (Ct) value higher than 30 were considered negative [40]. To confirm the initial putative findings of F. tularensis (Tul4+ samples), subsequent qPCRs were performed using the Francisella qPCR assays in Table 1 [40][41][42][43]. In brief, the DNA was pre-amplified (due to the limited amount of DNA) using the RepliG midi kit (Qiagen, Hilden, Germany), according to the manufacturer's instructions. The PCRs (25 µL) consisted of 2X PerfeCTa qPCR ToughMix (VWR, Radnor, PA, USA), 0.5 µM (final concentration) of each primer, 0.1 µM (final concentration) of each probe, and 1 µL template. The temperature profile was as follows: 95 • C for 10 min, followed by a two-step cycle of 15 s at 95 • C and 60 s at 60 • C. Positive and negative controls were included.
The tick DNA samples were tested for PCR inhibitor with an additional qPCR using seal herpesvirus type 1 [44].

Spotted Fever Group Rickettsia
Ticks also were screened for SFGR DNA by microfluidic qPCR, using primers and probes targeting the gltA gene, according to Michelet et al. [40]. Samples with a Ct-value higher than 30 were considered negative [40]. Confirmation analyses were executed on a small set of ticks (n = 38) with Ct-values gltA ranging from 5.6 to 29.6, using primers targeting the 17 kDa gene by Carl et al. [45] (Table 1). This was done because the confirmation PCR did not include a pre-amplification step, as was the case with the microfluidic qPCR, making it a less sensitive method. In brief, the reaction comprised of 1X Phusion Green HF buffer (

Tick Taxon
Tick taxa and Hyalomma species were determined by PCR, using primers by Beati and Keirans [43] (Table 1). See Hoffman et al. 2021 [39] for further details. Morphological determination was not performed since identification of species level of immature ticks belonging to the tick complex H. marginatum, including the species H. rufipes and H. marginatum, is difficult and not recommended [46]. Furthermore, life stage determination was not performed for the infesting ticks. However, the majority of the avian-associated ticks were likely immatures [36,47].

Sanger Sequencing
12S rDNA and 17 kDa amplicons were treated with illustra ExoProStar TM 1-step kit (Cytiva, Marlborough, MA, USA), according to the instructions by the manufacturer, prior Sanger sequencing at Macrogen (Amsterdam, the Netherlands).

Spotted Fever Group Rickettsia
The CLC Main Workbench 7 by Qiagen (Aarhus, Denmark) was used for assembling partial 17 kDa sequences, which were compared to sequences deposited in the GenBank database [48] using the nucleotide Basic Local Alignment Search Tool (BLASTN) (v.2.10.0) [49].

Metagenomic Sequencing
Two DNA samples putatively positive for F. tularensis by microfluidic qPCR were whole genome amplified (RepliG midi kit, Qiagen, Hilden, Germany) before being subjected to further characterization using two sequencing technologies. Non-enriched samples were prepared using TruSeq PCR-free library kits (Illumina, San Diego, CA, USA) and sequenced as 2 × 150 base pairs (bp) in one lane on an S4 flow cell in a NovaSeq 6000 sequencing instrument (Illumina) at the SNP&SEQ platform at NGI Uppsala (Sweden). Enriched samples were prepared using Nextera library kits (Illumina) and sequenced as 2 × 150 bp on a 300-cycle sequencing flow cell using a NextSeq instrument (Illumina) at NAU (Northern Arizona University). Non-enriched samples were also prepared using LSK-109 library kits (Oxford Nanopore Technologies, Oxford, UK) and sequenced in MinION R9.4.1 flow cells using a MinION sequencing device (Oxford Nanopore Technologies).

Enrichment
Due to few metagenomic reads mapping to Francisella, RNA baiting was performed at NAU to enrich Francisella DNA present in the tick extracts. Briefly, pre-amplified tick DNA extracts were uniquely indexed in~300 bp sequencing libraries and exposed to RNA hybridization baits (probes) of 120 bp (Agilent Technologies Inc., Santa Clara, CA, USA). The RNA baits were designed against a Francisella pan-genome defined by 498 Francisella genomes examined in [7]. Sequences <120 bp were removed, and regions with a homology of ≥80% with non-Francisella bacteria and ribosomal RNA genes were excluded, yielding 188,430 unique probe signatures, which included 2X tiling to ensure 50% sequence overlap to optimize capture. To further refine optimal capture, manufacturing of replicate copies of~20,000 probes comprised of high (≥50%) or low (≤22%) GC content was performed. Bait capture was executed twice for increased purification of Francisella sequences from background tick DNA.

Taxonomic Classification of Sequence Reads
Sequenced tick samples were characterized with a custom-made database containing bacteria, eukaryotes, and viruses, using Kraken 2 [50]. The database was created using FlexTaxD [51] with bacteria based on the taxonomy from the Genome Taxonomy Database (GTDB) together with eukaryotes and viruses from the National Center for Biotechnology Information (NCBI). The post-processing tool StringMeUp (v.0.1.4) (https://github.com/ danisven/StringMeUp, accessed on 15 May 2021) was used for adjusting the results for different confidence scores.

Tick Species Confirmation
Species confirmation of the metagenomic characterized ticks was performed using assembled mitochondrial genomes and the animal identification engine provided by BOLD [58], in which the mitochondrial cytochrome oxidase subunit 1 (COI) gene was used.

Phylogenetic Analyses Tick Phylogenies
Assembly of 12S rDNA sequences was performed in the CLC Main Workbench 7 (Qiagen, Aarhus, Denmark). Partial 12S rDNA sequences were aligned using the MAFFT algorithm and compared to sequences available in GenBank [48] using BLASTN (v.2.10.0) [49] and to sequences from morphologically determined reference specimens of multiple species of Hyalomma. Maximum likelihood 12S rDNA phylogenies were built in MEGA7 [59], and tick sequences were grouped based on their position in the 12S rDNA phylogenies. See Hoffman et al. 2021 [39] for details.

Biotin Synthesis Pathways
Previous analyses have identified multiple genes involved in the biotin synthesis pathway in the FLE of the tick O. moubata, including bioA, bioB, bioC, bioD, and bioF [11,13]. Homologs to these genes are present in the genome of the FLE of the tick species Argus arboreus (Francisella persica) [66]. To assess the conservation of these same genes in the genomes of the Francisella-like and Midichloria endosymbionts from H. rufipes samples D14IT15.2 and D14IT20, sequencing reads from both the enriched and non-enriched metagenomic data for these samples were mapped to these genes in the F. persica and M. mitochondrii genomes with minimap2 (v.2.22) [67] and the breadth of coverage was calculated with Samtools (v.1.11) [68] at a minimum depth of 3X.

Bird Trapping and Tick Collection
In total, 10,209 birds were trapped and screened for ticks. Of these, 244 (2.4%) birds were found to be infested by ticks (n = 575) ( Table S2 in the Supplementary Materials). Most of the tick-infested birds were long-distance migrants (98.0%). See Hoffman et al. [39] for additional details and information about the distribution pattern of ticks on the bird species.

Tick Determination
The collected ticks were assigned to the Ixodidae genera Hyalomma, Ixodes, Amblyomma, and Haemaphysalis, according to their position in phylogenies based on partial 12S rDNA sequences [39]. The assignment was not possible for 11.5% of the ticks due to the absence of PCR amplicons. The most common tick species were H. rufipes and H. marginatum [39]. The two metagenomics characterized ticks were confirmed as H. rufipes based on their similarity to a known COI sequence from H. rufipes and their positioning in the Hyalomma 12S rDNA [39] and mitochondrion (Figure 1) phylogenies.

Spotted Fever Group Rickettsia
The microfluidic qPCR data suggested the presence of SFGR in 59.1% (340/575, gltA+) of the total collected ticks, including 60.0% (

Metagenomics
Classification of metagenomic sequencing reads from tick samples that were putatively positive for F. tularensis by microfluidic qPCR confirmed the presence of FLEs and not F. tularensis. However, assembly of complete Francisella genomes was not possible from these data due to low read count and coverage (Table 3). Two rounds of Francisella enrichment were therefore performed, resulting in >95% Francisella DNA after enrichment. The assembly of Francisella DNA present in the sequenced enriched sample D14IT15.2 was in total 1,413,985 bp divided into 510 contigs with N50 = 10,189 and N50n = 36. The assembly of Francisella DNA present in the sequenced enriched sample D14IT20 was, in total, 1,415,455 bp divided into 506 contigs with N50 = 11,392 and N50n = 35.    Table 4. Highest average nucleotide identity for bacterial (n = 6) and tick (n = 2) metagenomeassembled genomes generated in this study detected in two avian-associated Hyalomma rufipes ticks (D14IT15.2 and D14IT20) that tested positive for Francisella tularensis and spotted fever group Rickettsia during screening.  1 According to the Genome Taxonomy Database; 2 According to the National Center for Biotechnology Information; 3 Mitochondrion; MAG-metagenome-assembled genome; ANIb-average nucleotide identity, BLAST method; FLE-Francisella-like endosymbiont; Om-Ornithodoros moubata.

Biotin Gene Conservation
Based upon the mapping of reads to the corresponding F. persica coding DNA sequences (CDSs) encoding homologs to bioA, bioB, bioC, bioD, and bioF, at least bioA appears to be missing in the FLE MAGs generated from the enrichment of samples D14IT15.2 and D14IT20, indicating the biotin pathway is not intact in these FLEs (Table 5). We note that mapping of FLE reads from the non-enriched metagenomic sequencing data to these same CDSs was very limited, owing to the low proportion of FLE reads in these data (Table 3). In contrast, there was significant mapping of reads from the non-enriched metagenomic data to all of the examined biotin genes in the M. mitochondrii genome, indicating that this pathway is likely intact in the Midichloria endosymbiont of H. rufipes. There was a limited mapping of reads from the enriched metagenomic data to the biotin genes in M. mitochondrii, which is not unexpected given the high proportion of Francisella DNA in these samples following enrichment.

Discussion
In this study, screening of ticks infesting birds migrating in a northward route from wintering areas in Africa revealed a high prevalence of Francisella in the tick species H. rufipes (fopA+: 76.7%), a known vector of SFGR and Crimean-Congo hemorrhagic fever virus [69], which suggest that migratory birds in the AWPR may contribute to northward dispersal of Francisella-infected ticks. The high prevalence of Francisella likely represents a high prevalence of FLEs, as FLEs have previously been detected in multiple species of Hyalomma ticks and at a high prevalence [17,20,70]. The two Francisella MAGs generated from two H. rufipes (Hr) were found to be members of a subclade within the FLE group (GTDB cluster: Francisella sp002095075) in Clade 1 of Francisella (Figure 2), which includes the FLE species FLE-Om (present in the soft tick O. moubata) and FLE-Am (present in the hard tick A. maculatum). The FLE-Hr had the highest identity to FLE-Om (ANIb: 96.7-97.0%). Several ticks tested positive for both Francisella and SFGR spp. (fopA+/gltA+: 47.1%), suggesting that presence of Francisella did not prevent the occurrence of SFGR spp.; a similar observation has been reported by Scoles [19]. ANI and phylogenetic inference indicated the highest similarity of the Rickettsia MAGs detected in H. rufipes (Rickettsia-Hr) with R. aeschlimannii (ANIb: 98.8-99.9%), a SFGR (i) reported to cause human infections [71,72], (ii) associated with Hyalomma ticks, including H. rufipes and H. marginatum [34,[73][74][75][76][77], (iii) previously identified together with FLEs in H. marginatum [23], and (iv) detected at similar prevalences as in this study in ticks of the H. marginatum species complex infesting northbound migratory birds [36,78]. Data from the two metagenomic sequenced H. rufipes ticks indicated co-occurrence also with a species of Midichloria. ANI and phylogenetic inference of the Midichloria MAGs detected in H. rufipes (Midichloria-Hr) indicated a close relationship to M. mitochondrii (ANIb: 91.5-92.3%). Midichloria sp. bacteria have previously been detected in H. marginatum species complex ticks infesting northbound trans-Saharan spring migrating birds trapped in Italy [47]. That study found a high prevalence of Midichloria DNA in the investigated Hyalomma ticks (>90%) and in a considerable fraction of the blood samples from the avian hosts (>40%) and suggested that the presence of Midichloria DNA in the blood was associated with lower fat reserves in the tick-infested birds [47].
Ticks may depend on FLEs because they provide nutrients that are absent in the tick diet, such as B vitamins (folate/folic acid (B 9 ), riboflavin (B 2 ), and biotin (B 7 )) and co-factors, and thereby improve the fitness of the tick [5,11]. It has been suggested that FLEs serve as alternative obligate symbionts in some species of ticks [5,79], whereas Coxiella-like endosymbionts (CLEs) are considered to be obligate symbionts (i.e., present in most specimens) in most tick species [6,80]. Ticks may escape a negative symbiosis by replacing an old symbiont with a new bacterium [81]. FLEs may have replaced CLEs in several tick lineages, including O. moubata and A. maculatum [5,13,79]. M. mitochondrii has also been suggested to be a nutritional endosymbiont since it encodes genes for the production of several co-factors and B vitamin biotin [82]. Buysse et al. [23] showed that the genomes of FLEs detected in H. marginatum included functional biosynthesis pathways for folate and riboflavin but were deprived of a functional biosynthesis pathway for biotin. The authors suggested that this was compensated for by the co-symbiosis with Midichloria bacteria also present in H. marginatum, since their genomes included an intact biotin biosynthesis operon [23]. The Midichloria detected in H. marginatum had a partial riboflavin biosynthesis pathway, indicating that co-occurrence of FLEs and Midichloria may be essential for complete nutritional symbiosis in H. marginatum [23]. We observed similar patterns for H. rufipes: a disrupted biotin biosynthesis pathway in its FLE but an apparently intact biotin biosynthesis pathway within its Midichloria endosymbiont (Table 5), suggesting the co-occurrence of these bacteria also may be essential for nutritional symbiosis in H. rufipes. As noted by Buysse et al. [23], a similar dual symbiosis may be present in several other Hyalomma tick species (Hyalomma aegyptium, Hyalomma anatolicum, Hyalomma dromedarii, Hyalomma excavatum, Hyalomma impeltatum, Hyalomma lusitanicum, and Hyalomma truncatum), as they found evidence for the presence of both FLEs and Midichloria within them. The apparent exception within this tick genus to date appears to be Hyalomma asiaticum, which does not harbor any Midichloria but instead has an FLE with an intact biotin pathway [83].
Molecular species determination of members of the H. marginatum species complexcurrently consisting of five species, including H. marginatum and H. rufipes [84]-can be difficult due to the inclusion in public databases of sequences obtained from incorrectly identified specimens [46]. Complete tick mitochondrial MAGs were therefore constructed to verify the initial 12S rDNA-based speciation of the metagenomically characterized ticks. The two characterized tick specimens were found to group in the H. rufipes clade also in the mitochondrion phylogeny, verifying the initial 12S rDNA speciation results.

Conclusions
Understanding the biology, ecology, and evolution of tick endosymbionts is important, as they may share a close evolutionary relationship with pathogenic bacteria and also may influence the fitness [4] and even the behavior of the tick host [85]. The results of this study demonstrate that FLEs are present in many H. rufipes ticks, and migratory birds in the AWPR contribute to the northward geographical spread of FLE-containing ticks. The absence of F. tularensis in the investigated ticks does not provide evidence supporting that immature life stages of H. rufipes contribute to the transmission of F. tularensis in the study region. Furthermore, the results suggest that migratory birds also contribute to northward geographical spread in the AWPR of H. rufipes ticks containing SFGR spp., including R. aeschlimannii and Midichloria bacteria, and that a dual endosymbiosis (co-symbiosis) of FLEs and Midichloria may support the nutritional requirements of the medically important tick vector H. rufipes. We acknowledge that the majority of the results of this study are based on unconfirmed screening data and that the reported detection results are therefore conservative estimates of the prevalence. Future studies should therefore focus on verifying the Francisella and SFGR prevalence in H. rufipes as well as investigate the Midichloria prevalence in H. rufipes and the impact that FLEs and Midichloria may have on H. rufipes, including their interaction with bacterial pathogens, such as SFGR.
Supplementary Materials: Supporting information can be downloaded at: https://www.mdpi. com/article/10.3390/microorganisms10071393/s1, Table S1. Publicly available genomes for Francisella, Rickettsia, Midichloria, and tick mitochondria downloaded from GTDB and NCBI. Table S2. Microfluidic real-time PCR data for ticks collected from bird species trapped in the Mediterranean basin during the spring migration of 2014 and 2015. The ticks were screened for Francisella and spotted fever group Rickettsia species using primers targeting the fopA and gltA genes. Figure S1. Whole genome maximum likelihood phylogeny. Highlighted area indicates Midichloria clade. Study genomes are in bold. Orientia tsutsugamushi strain Karp was used to root the tree. Bootstrap values ≥ 75 are presented at the nodes. The scale bar represents the expected number of substitutions per site. Figure S2. Whole genome maximum likelihood phylogeny of Rickettsia. Highlighted area indicates clade for Rickettsia aeschlimannii. Study genomes are in bold. Orientia tsutsugamushi was used to root the tree. Bootstrap values ≥ 75 are presented at the nodes. The scale bar represents the expected number of substitutions per site. Figure S3. Heatmap of the average nucleotide identity (ANI), demonstrating nucleotide-level genomic similarity between Francisella genomes. The pairwise comparison of 19 Francisella genomes was computed by BLAST, using the pyANI software. Study genomes are in bold. FLE, Francisella-like endosymbiont. Figure S4. Heatmap of the average nucleotide identity (ANI), demonstrating nucleotide-level genomic similarity between Rickettsia genomes. The pairwise comparison of 138 Rickettsia genomes was computed by BLAST, using the pyANI software. Study genomes are in bold. For NCBI organism names, see Table S1. Figure S5. Heatmap of the average nucleotide identity (ANI), demonstrating nucleotide-level genomic similarity between bacterial genomes. The pairwise comparison of 23 genomes was computed by BLAST, using the pyANI software. Study genomes are in bold. Figure S6. Heatmap of the average nucleotide identity (ANI), demonstrating nucleotide-level genomic similarity between Hyalomma mitochondrial genomes. The pairwise comparison of 23 Hyalomma mitochondrial genomes was computed by BLAST, using the pyANI software. Study genomes are in bold.

Data Availability Statement:
The metagenomic sequence data generated in this study are available under the NCBI BioProjectPRJNA764565.