# Deep sequencing of 16 Ixodes ricinus ticks unveils insights into their interactions with endosymbionts

Paulina Lesiczka, Tal Azagi, Aleksandra Krawczyk, William Scott, Ron Dirks, Ladislav Simo, Gerhard Dobler, Bart Nijsse, Peter Schaap, Hein Sprong, et al.

► To cite this version:

Paulina Lesiczka, Tal Azagi, Aleksandra Krawczyk, William Scott, Ron Dirks, et al. Deep sequencing of 16 Ixodes ricinus ticks unveils insights into their interactions with endosymbionts. mSystems, 2025, pp.e00507-25. 10.1128/msystems.00507-25. hal-05132420

HAL Id: hal-05132420

https://hal.inrae.fr/hal-05132420v1

Submitted on 21 Jul 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Deep sequencing of 16 *Ixodes ricinus* ticks unveils insights into their interactions with endosymbionts

Paulina M. Lesiczka,¹,² Tal Azagi,³ Aleksandra I. Krawczyk,³ William T. Scott, Jr.,⁴,⁵ Ron P. Dirks,⁶ Ladislav Simo,⁷ Gerhard Dobler,⁸ Bart Nijsse,⁴,⁵ Peter J. Schaap,⁴,⁵ Hein Sprong,³ Jasper J. Koehorst⁴,⁵

**AUTHOR AFFILIATIONS** See affiliation list on p. 17.

**ABSTRACT** *Ixodes ricinus* ticks act as vectors for numerous pathogens that present substantial health threats. Additionally, they harbor vertically transmitted symbionts, some of which have been linked to diseases. The difficulty of isolating and cultivating these symbionts has hampered our understanding of their biological role, their potential to cause disease, and their modes of transmission. To expand our understanding of the tick symbiont *Midichloria mitochondrii* and *Rickettsia helvetica*, which has been linked to disease in humans, we utilized deep sequencing on 16 individual adult female ticks collected from coastal dune and forested areas in the Netherlands. By employing a combination of second- and third-generation sequencing techniques, we successfully reconstructed the complete genomes of *M. mitochondrii* from 11 individuals, *R. helvetica* from eight individuals, and the mitochondrial genome from all ticks. Additionally, we visualized the location of *R. helvetica* in tick organs and constructed genome-scale metabolic models (GEMs) of both symbionts to study their environmental dependencies. Our analysis revealed a strong cophylogeny between *M. mitochondrii* and mitochondrial genomes, suggesting frequent maternal transmission. In contrast, the absence of cophylogeny between *R. helvetica* and the mitochondrial genomes, coupled with its presence in the receptaculum seminis of *I. ricinus* females, raises the possibility of paternal transmission of *R. helvetica*. Notably, the genetic diversity of *R. helvetica* was found to be very low, except for the rickA virulence gene, where the presence of up to 13 insertions of a 33 nt-long repeat led to significant variability. However, this variation could not account for the differences in infection prevalence observed across eight distinct locations in the Netherlands. By employing deep sequencing, it becomes feasible to extract complete genomes and genetic data of symbionts directly from their host organisms. This methodology serves as a robust means to gain fresh insights into their interactions. Our observations, which suggest paternal transmission of *R. helvetica*, a relatively unexplored mode of transmission in ticks, require validation through experimental investigations. The genetic variations identified in the rickA virulence gene of *R. helvetica* have the potential to influence the infectivity and transmission dynamics of *R. helvetica*.

**IMPORTANCE** Ticks are vectors of numerous human pathogens; however, the microbial interactions within ticks and the mechanisms governing pathogen transmission remain poorly understood. This study uses deep sequencing of individual *Ixodes ricinus* to reconstruct high-quality genomes of endosymbionts and the mitochondrion of the tick, revealing previously undetected microbial dynamics. Notably, we recovered low-abundance *Rickettsia* and *Midichloria* genomes from single ticks and present evidence that suggests paternal transmission of *R. helvetica*. These findings offer novel insights into the ecology and evolution of tick-associated microbes and have implications for understanding the origins and spread of tick-borne diseases.

Editor Jonathan L. Klassen, University of Connecticut, Storrs, Connecticut, USA

Address correspondence to Hein Sprong, hein.sprong@rivm.nl, or Jasper J. Koehorst, jasper.koehorst@wur.nl.

Paulina M. Lesiczka and Tal Azagi contributed equally to this article. The author order was determined by drawing straws.

The authors declare no conflict of interest.

See the funding table on p. 17.

Received 11 April 2025

Accepted 16 May 2025

Published 16 June 2025

Copyright © 2025 Lesiczka et al. This is an open-access article distributed under the terms of the Creative Commons Attribution 4.0 International license.
**KEYWORDS** *Ixodes ricinus*, deep sequencing, symbionts, genome reconstruction, *Midichloria mitochondrii*, *Rickettsia helvetica*, genome-scale metabolic modeling, cophylogenetic analysis, paternal transmission

*Ixodes ricinus*, the most common tick species in the Palearctic realm, carries a variety of facultative symbionts, some of which can cause infectious diseases in humans. The most common tick-borne diseases associated with *I. ricinus* are Lyme disease (LB) (which is caused by *Borrelia burgdorferi sensu lato*) and tick-borne encephalitis (TBE), which is caused by tick-borne encephalitis virus (1). In addition, *I. ricinus* transmits several emerging pathogens such as *Rickettsia helvetica*, *Anaplasma phagocytophilum*, *Borrelia miyamotoi*, *Neoehrlichia mikurensis*, and several *Babesia* species (2). Infections with these microorganisms have been associated with human diseases in cross-sectional and case studies, but their pathogenic potential has not been fully understood (3, 4). These microorganisms can influence the tick host's metabolism, development, reproduction, stress defense, and immunity (5) and may indirectly contribute to the transmission dynamics of tick-borne pathogens, significantly impacting human and animal health (6). This implies a broader role within the tick host, and thus, we refer to them as tick symbionts.

Bacterial endosymbionts of ticks are bacteria living within the body or cells of ticks, and their survival is critically dependent on the life cycle of ticks. Some of those bacteria might also have the ability to cause infection in vertebrate hosts, more specifically humans, thereby causing disease. Such bacteria can therefore also be considered (human) pathogens. In other words, "pathogen" and "symbiont" are not mutually exclusive properties/descriptions of bacteria living inside ticks (7).

The biological role, virulence, and transmission dynamics of tick symbionts and some of the pathogens are often difficult to study due to challenges in isolating and culturing them. A notable example is *Midichloria mitochondrii*, initially discovered in the ovaries of ticks. *Midichloria mitochondrii* has been found to be highly prevalent in various tissues of *I. ricinus* (8–11). Since its discovery, *Midichloria* symbionts have been detected in several tick species across different genera, albeit with varying prevalence (12), and are believed to contribute to detoxification, vitamin supply, and energy-related functions, based on their observed interactions and effects within host cells (13). Although *M. mitochondrii* has been detected in the blood of animals, it has yet to be associated with human disease and is mainly considered in the context of tick-symbiont interactions (14), and other symbionts such as *Rickettsia helvetica* cannot be so clearly defined.

*R. helvetica* is a gram-negative, obligate intracellular bacterium formerly included in the spotted fever group rickettsiae (SFGR) (15). *R. helvetica* can be transmitted by tick bites and has been associated with human disease. Only a few human cases have been reported in Europe (16–21), and detection in cerebrospinal fluid of suspected Lyme neuroborreliosis patients has been shown in the Netherlands and Denmark. However, it is unclear whether infection with *R. helvetica* can cause disease or might change the course of Lyme disease (22, 23). Moreover, studies on the involvement of *R. helvetica* in sarcoidosis seem to be contradictory (24). Thus, the extent of *R. helvetica*'s infectivity and pathogenicity is still under investigation, and little is known about its ecoepidemiology (25, 26). DNA of *R. helvetica* has only been sporadically detected in wildlife (27–31) with no clear animal reservoir. This symbiont is transmitted vertically from one generation of ticks to the next and has been repeatedly isolated from questing *I. ricinus* in European countries (32), which are currently thought to be the main reservoir for *R. helvetica* in nature (27). *R. helvetica* infection rates in *I. ricinus* vary from 0% to 66%, resulting in highly variable human exposure from tick bites. The reason for this high variability is unknown and could be attributed to the capacity of *R. helvetica* to be beneficial for the tick host under certain environmental conditions (33), or it could be related to genetic variation within *R. helvetica* populations, but the latter has not been investigated so far.

The pathogenic potential of *R. helvetica* has been explored, to a certain extent, on the cellular level. Although other rickettsial species, which are pathogenic to humans, utilize
actin-based cell-to-cell motility (34), a culture of *R. helvetica* has been shown to lack this ability *in vitro* (26), possibly due to a deletion in the *sca2* gene, encoding a protein required for virulence that stimulates actin filament assembly (35, 36). In addition, although the sequence encoding the RickA protein, a bacterial actin nucleator (37), was complete, its expression was apparently insufficient to enable actin-based motility in the strains studied (26). Comparative analyses of its genetic makeup and metabolic potential are needed to clarify the pathogenic potential of *R. helvetica* as well as its interaction and localization within the tick host.

Advances in sequencing technologies have ushered in a new era of understanding by enabling the comprehensive annotation of complete tick mitochondrial genomes (mt genomes). These genomic insights have significantly enriched our understanding of tick phylogeny, population genetics, and evolutionary processes (38, 39). Given ticks' intimate association with symbiotic microbes crucial for their hematophagous lifestyle (40), the incorporation of mt genomes into co-phylogeny analyses presents a promising avenue for unraveling the evolutionary patterns of tick symbionts (41). This approach holds particular significance for symbionts like *M. mitochondrii*, which inhabit tick mitochondria offering a unique window into their evolutionary trajectory (42).

Genome-scale model (GEM) reconstruction and analysis can further enhance this understanding by detailing the metabolic interactions and dependencies between ticks and their symbionts. These models can elucidate how these endosymbionts may contribute to the metabolic network, revealing potential targets for controlling tick populations and preventing tick-borne diseases. Additionally, GEMs can provide insights into the co-evolutionary dynamics of ticks and their symbionts, as has been shown in recent work in insects (43), offering a clearer picture of how these complex biological systems have adapted and evolved together over time.

In this study, we applied deep sequencing on 16 *l. ricinus* specimens to unveil the genetic variation and differences in transmission dynamics of *M. mitochondrii* and *R. helvetica*. Complete genomes were obtained directly from the individual ticks using a hybrid assembly approach combining PromethION and Illumina sequencing. In addition to interspecies whole genome comparisons, we explored the genetic association of *R. helvetica* and *M. mitochondrii* with their tick hosts to better understand the mode of transmission of these symbionts. To further our understanding of the latter, we visualized *R. helvetica* in tick organs.

# MATERIALS AND METHODS

## Tick specimen metadata, DNA isolation and sequencing, and data FAIRification

### Tick specimen metadata and processing

#### Field samples

Six questing *l. ricinus* adult females, Ir_f1, Ir_f3, Ir_f6, Ir_f11, Ir_f12, and Ir_f16, were collected by blanket dragging in a planted forest region (f) in the Netherlands. Another eight adult females were collected in a coastal sand dune coastal area (d) in the Netherlands and labeled as Ir_d1, Ir_d2, Ir_d3, Ir_d4, Ir_d5, Ir_d6, Ir_d7, and Ir_d9. The exact locations are in Table 1. Ticks were kept at -80°C until further processing.

#### Lab-reared adult females

Questing *l. ricinus* nymphs were collected by blanket dragging and morphologically identified to species level using an identification key (44). Nymphs were fed on an artificial membrane blood-feeding system, as described in reference 45. Two lab-reared adult females, F1 and F2, were used in this study (Table 1).
<table><caption>TABLE 1 Origin and deep sequencing of adult *I. ricinus* ticks</caption><thead><tr><th>Description</th><th>Collection date</th><th>Location (habitat type)</th><th>Name</th><th>ONT (Gb)</th><th>Illumina (Gb)</th></tr></thead><tbody><tr><td rowspan="8">Questing adult female</td><td rowspan="8">14 October 2021</td><td rowspan="8">52.3037; 4.5331 (coastal dune)</td><td>Ir_d1</td><td>65.7</td><td>157.0</td></tr><tr><td>Ir_d2</td><td>19.4</td><td>172.4</td></tr><tr><td>Ir_d3</td><td>16.4</td><td>91.4</td></tr><tr><td>Ir_d4</td><td>14.4</td><td>156.2</td></tr><tr><td>Ir_d5</td><td>14.9</td><td>115.8</td></tr><tr><td>Ir_d6</td><td>32.9</td><td>192.0</td></tr><tr><td>Ir_d7</td><td>48.3</td><td>169.9</td></tr><tr><td>Ir_d9</td><td>38.4</td><td>158.5</td></tr><tr><td rowspan="7">Questing adult female</td><td rowspan="7">23 March 2022</td><td rowspan="7">51.9769; 5.6992 (planted forest)</td><td>Ir_f1</td><td>144.1</td><td>137.7</td></tr><tr><td>Ir_f3</td><td>88.4</td><td>105.5</td></tr><tr><td>Ir_f6</td><td>69.7</td><td>111.3</td></tr><tr><td>Ir_f11</td><td>63.7</td><td>140.6</td></tr><tr><td>Ir_f12</td><td>39.0</td><td>120.0</td></tr><tr><td>Ir_f16</td><td>67.1</td><td>139.1</td></tr><tr><td>Questing nymph lab-reared to female</td><td>9 January 2020</td><td>51.5839; 5.4159 (planted forest)</td><td>F1</td><td>53.0</td><td>129.8</td></tr><tr><td>R. helvetica culture from questing male</td><td>14 August 2018</td><td>52.4347; 4.6258 (coastal dune)</td><td>F2</td><td>94.8</td><td>104.8</td></tr><tr><td>R. helvetica culture from questing nymph</td><td>2 July 2013</td><td>48.3177; 12.5390 (planted forest)</td><td>DK2</td><td>1.22</td><td>0.63</td></tr><tr><td></td><td></td><td></td><td>OB144/2013</td><td>8.71</td><td>4.55</td></tr></tbody></table>

## DNA isolation and sequencing of biospecimens

DNA extraction and sequencing were performed at Future Genomics Technologies BV (Leiden, The Netherlands). Ticks were ground to a fine powder using a pestle and mortar with liquid nitrogen. To extract high molecular weight DNA, a Genomic-tip 20 /G kit (Qiagen Benelux BV, Venlo, The Netherlands) was used as per protocol.

DNA quality was measured via electrophoresis in Genomic DNA ScreenTape on an Agilent 4200 TapeStation System (Agilent Technologies Netherlands BV, Amstelveen, The Netherlands) and total DNA using a Qubit 3.0 Fluorometer (Life Technologies Europe BV, Bleiswijk, The Netherlands). Prior to sequencing, detection of *R. helvetica* was done by qPCR of the *gltA* gene using *Rick_HelvgltA_F2*: ATGATCCGTTTAGGTTAATAGGCTTCGGTC and *Rick_HelvgltA_R2*: TTGTAAGAGCGGATTGTTTTCTAGCTGTC as forward and reverse primers as described (46) (File S1.1).

One microgram of the DNA samples was used to prepare a 1D ligation library using the Ligation Sequencing Kit (SQK-LSK110) as per protocol (Oxford Nanopore Technologies [ONT], Oxford, United Kingdom). ONT libraries were tested on a MinION flowcell (FLO-MIN106) and subsequently run on a PromethION flowcell (FLO-PRO002) using the following settings: basecall model: high-accuracy; basecaller version: Guppy 5.0.17. Parallel aliquots of the DNA samples used for ONT sequencing were used to prepare Illumina libraries using the Nextera DNA Flex Library Prep Kit as per protocol (Illumina Inc., San Diego, CA, USA). Library quality was measured using an Agilent 4200 TapeStation System as described above. Genomic paired-end (PE) 150 nt libraries were sequenced using the Illumina NovaSeq 6000.

## DNA isolation and sequencing of *R. helvetica* strains from Vero cells

*R. helvetica* strain OB144 originated from a questing nymph in a planted forest in Germany, and DK2 from a questing adult male tick (Table 1). Strains were cultured and harvested from VERO cells as described (47, 48) with two modifications: here, the needle and syringe protocol was implemented using 26-gauge needles and a 0.45-µm syringe-driven membrane filter. Genomic DNA was extracted using the ZymoBIOMICS
96 MagBead DNA Kit (Zymo Research, Orange, CA). OB144 was sequenced by Future Genomics Technologies BV as described above. DK2 was sequenced by Baseclear (Leiden, the Netherlands) using a NovaSeq 6000 system (Illumina, San Diego, CA, USA). DK2 long reads were generated by RIVM using a Nanopore GridION sequencer using the sequencing kit LSK110. Base calling was performed using the ONT Guppy barcoding software version 4.4.1+1c81d62.

## Taxonomic classification of Illumina reads

Illumina reads were filtered with fastp v0.2.23 with parameters -correction, -disable_TRIM_POLY_g, -length_required 50, -qualified_quality_phred 20, -unqualified_percent_limit 20. For taxonomic classification of the Illumina reads, the filtered reads were run through Kraken2 v2.1.2 in paired-end mode (-paired-end) with a confidence threshold (-confidence) of 0.5. For the Kraken2 index, the full NCBI nucleotide (nt) database (downloaded from https://genome-idx.s3.amazonaws.com/kraken/k2_nt_20230502.tar.gz) was used. The index was constructed using default parameters: kmer length of 35, minimizer length of 31, and minimizer spacing of 7. For species abundance estimation, Bracken v2.8 was used with the previously described NCBI nt kraken2 database with a read length (-r) of 150 and a threshold (-t) of 1,000 reads and level (-l) "S" for species. Eukaryotic species were filtered from the Bracken output files with the "filter bracken out.py" script from KrakenTools v1.2e. To get read counts below 1,000 for the resulting species, Bracken was run on the same parameters except with a threshold (-t) of 10 reads.

## Genome assemblies and annotation

A standardized hybrid workflow for metagenomics assembly (49) was used, which is registered at the WorkflowHub (50). Illumina and nanopore reads were inspected for quality using FastQC v0.12.1 (51). Illumina reads were filtered with fastp v0.23.2 (52) with parameters -correction, -dedup, -disable_TRIM_POLY_g -length_required 50, -qualified_quality_phred 20, -unqualified_percent_limit 20. Possible Illumina spike in PhiX genomic reads contamination was filtered with BBMap/bbduk v38.95 (53) -k 31, reference (NC_001422.1). To filter Illumina reads for selected species, BBMap v38.95 (53) was used with default settings and bloom = t. Nanopore reads were filtered using minimap2 v2.24 (54) default settings, -x and map-ont. Mapped reads were kept. Genome sequences used for filtering were extracted from the European Nucleotide Archive. For *Rickettsia*, all available *Rickettsia* reference genomes were used. For *M. mitochondrii*, the reference genome was CP002130.1. For mitochondrial reads, all available *Ixodes* mitochondrial genomes were used. A list of genome accessions for filtering can be found in File S1.2.

## R. helvetica: rickA gene variability analyses

Questing nymphs were collected in Duin & Kruidberg, Schiermonnikoog, Gieten, Veldhoven, Montferland, and Wassenaar, the Netherlands as described (55), immersed in 70% alcohol and stored at -20°C. DNA was extracted by alkaline lysis as described (56). *R. helvetica* was detected by qPCR targeting the *gltA* as described above. In the next step, the *rickA* gene was detected by PCR assay targeting a 600 nt fragment of the gene on qPCR-positive samples, as described (8). DNA-free water was used as a negative control. Primers used were RIC1200 F2 GGCAAAATGTTAAAAATGTTT and RIC + 430 R3 CCRGYTTTTTAACCGTAGTAG. The *rickA* gene PCR products and the sequences from the new *R. helvetica* genomes were aligned with ClustalW using Geneious 11.1.4 software (57, 58).

## Phylogenetic analysis of *Ixodes ricinus* mitogenomic data

ClustalW alignments were computed using Geneious 11.1.4 software (57, 58). Phylogeny was inferred by IQTREE 2.1.3 (59), and the best-fit evolution model was selected based
on the Bayesian information criterion (BIC) computed by implemented ModelFinder (60). Branch supports were assessed by the ultrafast bootstrap (UFBoot) approximation (61) and the Shimodaira–Hasegawa-like approximate likelihood ratio test (SH-aLRT) (62). Tree was rooted at the midpoint. No out-group was included. The tree was visualized and edited in FigTree (v1.4.1) and Inkscape (v 0.91) (63).

## Cophylogenetic analysis

The single-copy core genes of 11 *R. helvetica* and 12 *M. mitochondrii* genomes were identified in Anvi'o 7 and made up of 1,484 and 1,127 single-copy gene (SCG) clusters, respectively. Gene alignments from SCG were extracted and concatenated using the program anvi-get-sequences-for-gene-clusters, with parameters: -concatenategene-cluster; -report-DNA-sequences. Maximum likelihood phylogenetic trees were constructed in IQ-TREE v2.2.0 (59). The best model according to the Bayesian information criterion (BIC) was selected with Model Finder (60) as implemented in IQ-TREE. All models were run with ultrafast bootstrapping (61). eMPRess, a program for phylogenetic tree reconciliation under the duplication-transfer-loss (DTL) model (64), was used to determine whether the mitochondrial genome of *I. ricinus* and *I. ricinus* symbionts: *M. mitochondrii* and *R. helvetica* coevolved by comparing their phylogenetic relationships. The input trees were based on the results of phylogenetic analyses. The analysis was conducted using the default eMPRess parameters.

## Genome-scale metabolic model reconstruction and analysis

Genome-scale metabolic models (GEMs) were developed using CarveMe (v1.5.1) (65), starting with their annotated genome sequences (Ir_d9Mm and Ir_d9Rh). For reaction essentiality analysis, we performed gene or reaction deletion analysis using flux balance analysis (FBA) (File S2). MEMOTE (v0.13.0) (66) was used to assess the quality of the GEMs. In addition, a full reaction essentiality analysis was performed on both GEMs, identifying all metabolites for which the symbionts rely on the host (File S2). Furthermore, we performed a permutation test to assess the significance of the overlap in essential reactions between the symbionts. COBRApy (v0.26.3) (67), employing Python (v3.11.3), was used for modeling simulation and analysis of each GEM using the IBM CPLEX Solver. Scripts used to reproduce the model simulations are provided on GitLab (https://git.wur.nl/unlock/projects/ticks).

## Immunofluorescence localization studies of *R. helvetica* in ticks

Fifteen *I. ricinus* questing female ticks were collected in the Amsterdamse Waterleidingduinen, the Netherlands, and morphologically identified to species level using an identification key (44). All ticks were kept alive until further processing. Ticks were dissected, and the gut, Malpighian tubules, salivary glands, reproductive system, and synganglion were collected. Tick carcasses were screened for the presence of *R. helvetica* DNA by qPCR. The same was performed on 10 control ticks from a laboratory colony (BIPAR, France).

Organs were fixated in 4% paraformaldehyde (PFA) in PBS for 2 h and then replaced with 70% ethanol and stored at -20C until used. The organs were washed with 0.5% Triton X-100 in PBS (PBST) five times, then replaced by 5% normal goat serum (Sigma) in PBST and incubated for 1 h at room temperature. Afterward, the serum was replaced by *R. conorii* commercial antibody 1:30 in PBST and incubated at 4°C for 2 days. Samples were incubated with DAPI before mounting in Prolong Antifade Diamond Mountant (Life Technologies), and images were obtained using a z-stack tool in Leica DMI8 confocal microscopy. Assembled z-stacks were adjusted in Adobe Photoshop version 24.0.1 (Adobe Inc.).
## Statistical analysis

Haplotype diversity (Hd), variance of Hd, nucleotide diversity ($\pi$), variable (polymorphic) sites, and parsimony informative sites were calculated for mitochondrial genomes of *I. ricinus* using the DNA SP6 (v6.12.03) (68).

To test whether the composition of repeats differs per prevalence and per area type (dune vs. forest), we performed PERMANOVA based on Bray-Curtis dissimilarity matrix using the adonis2 function from the vegan package (69).

## RESULTS

### Deep sequencing of *I. ricinus* adult females reveals the presence and diversity of symbiotic bacteria

Sixteen adult female specimens were sequenced using both Illumina and ONT (Oxford Nanopore Technologies) methodologies: eight from a coastal sand dune area (d), six from a planted forest area, and two specimens molted under laboratory conditions. Using Illumina 150 × 150 nt paired-end reads as input and the nt database (70) as a reference, we applied a Kraken2-Bracken workflow (46, 71) to estimate for each specimen the abundances and taxonomic classifications of the microbiome-related Illumina reads obtained (Fig. 1, step 1). In total, more than 90% of the read pairs that were attributed to the bacterial domain could be classified at the species level.

A significant number of paired-end reads (with cumulative scores of >20,000 read pairs) were assigned to known tick-borne pathogens: *Anaplasma phagocytophilum* was present in tick Ir_d1 and Ir_d9; *Borrelia afzelii* was present in tick Ir_d9; and *Neoehrli-chia mikurensis* (73) in tick Ir_f3. *Rickettsiella* species have recently been identified as facultative endosymbionts of *I. ricinus* ticks (74), and with more than 1 million read pairs assigned, *Rickettsiella* species appear to be present in Ir_d1, Ir_d6, and Ir_d7 (Fig. 2). *Rickettsia* was detected in eight samples, and in almost all cases, *R. asiatica* was the most prominent classified species, followed by *R. helvetica* and *R. conorii* (File S1.3).

### Genome reconstruction of *M. mitochondrii* endosymbiotic strains

Kraken analysis revealed that 15 ticks are positive for *M. mitochondria* (File S1.3). The exception was F1, raised in a laboratory setting (Fig. 2). High-quality MAGs of *M. mitochondrii* could be obtained from 11 ticks and assembly drafts from the other four (File S1.4).

The MAGs shared over 99% ANI with the *M. mitochondrii* reference strain (75) and were characterized by a high BUSCO score (>92.9%). Genome sizes ranged from 1,190,312 to 1,191,611 bp except for the assembly obtained from tick Ir_d2, which was 1,268,036 bp due to a duplication event (File S1.3 and S1.4). Overall coverage was between 24.6 and 103.10 (File S1.10). The average guanine-cytosine (GC) content

<figure><img src="image_2.png" /><figcaption>**FIG 1** *I. ricinus* endosymbiont sequence classification, genome assembly, and analysis. Illumina reads were used as input for the Kraken2-Bracken workflow (57, 72) to obtain microbiome-related abundances and taxonomic classifications. Subsequently, a standardized workflow for hybrid metagenomics assembly was used to obtain *Rickettsia* and *Midichloria* genomes from each positive specimen (see Materials and Methods for details). Two sources of genomic DNA were used as input: *I. ricinus* metagenomic DNA obtained from adult female ticks and *R. helvetica* genomic DNA isolated from infected Vero cell lines.</figcaption></figure>
<figure><img src="image_3.png" /><figcaption>FIG 2 Kraken2-Bracken distribution of classified Illumina paired-end reads at the species level. Default parameters for kmer length (35), minimizer length (31), and minimizer spacing (7) were used in this analysis. Left: top-level classification, right: classification of the bacterial domain at species level. The top three most prevalent species are indicated. A more detailed analysis is available in File S1.3.</figcaption></figure>

was 36.6%. All *M. mitochondrii* genomes contained three rRNAs, between 35 and 39 tRNAs, and the number of protein-coding genes varied from 1,322 to 1,448. The overall coverage of these regions was analyzed, and no deviations were observed, indicating that an assembly error is unlikely.

Among the eleven *de novo* *Midichloria* genome assemblies, a high level of genomic collinearity was observed (Fig. 3, panel B, File S3); however, in comparison with the reference strain, a significant amount of structural variation was observed. Genome variant calling revealed strain variations. The total number of variants ranged from 604 to 964 (8,509 total). Of these, 5,029 single-nucleotide polymorphisms (SNPs), 2,246 insertions, 1,989 deletions, and 335 complex variants were detected (Table 2). Although these variants likely reflect biological differences, the potential influence of (454) sequencing artifacts on insertion and deletion counts cannot be fully excluded from the interpretation.

EggNOG-based functional analysis suggests that the reference genome, and the 11 newly sequenced genomes encode a (near) complete folate biosynthetic pathway (six of the seven required genes detected) and a likely complete biotin biosynthetic pathway (all required genes detected). COG0331 is present in all newly sequenced genomes but absent from the reference genome. COG0331 corresponds to malonyl CoA-acyl carrier protein transacylase, encoded by the *fadD* gene, and is responsible for initiating the type II fatty acid synthesis pathway (76). A complete eggNOG-based functional analysis is provided in File S1.5.

## Genome reconstruction of *R. helvetica* endosymbiotic strains

Of the 14 specimens collected in the field (Table 1), eight tested positive for the presence of *R. helvetica* DNA by qPCR analysis (File S1.1). Kraken analysis using the complete NT-database as a reference, however, suggested that *R. asiatica*, *R. conorii*, and *R. helvetica*, three species that share almost identical genome sequences, appeared to exhibit a significant concurrent presence in these ticks (File S1.9). We used FastANI (77) to gain more insight into their evolutionary distances. A commonly used operational definition for prokaryotes states that if two genomes have an average nucleotide identity (ANI) of over 95%, they belong to the same species. Between them, the reference strains of *R. helvetica* and *R. asiatica* show >98.66% ANI and in cross-comparison with *R. conorii*, >93.3% ANI. Although these high identity scores explain the Kraken
<figure><img src="image_4.png" /><figcaption>FIG 3 Major intragenomic rearrangements between *R. helvetica* and *R. asiatica*. FastANI orthologous mappings were done using a 3 kb fragment size. Red line segments denote the orthologous mappings between these fragments. (A) ANI computed between *R. conorii* strain (NC_003103.1) as query and *R. helvetica* strain C9P9 (NZ_CM001467.1) as reference. (B) ANI computed between *R. asiatica* strain (NZ_AP019563) as query and *R. helvetica* strain C9P9 as reference. (C) ANI computed between the *de novo* assembled *Rickettsia* strain from tick *Ir_dr9* as query and *R. helvetica* strain C9P9 as reference.</figcaption></figure>

results, it is important to note that FastANI revealed significant structural variation (Fig. 3). An ANI analysis between the reference strain *R. helvetica* and the reconstructed genomes showed an ANI value >99.9% and >98.6% for *R. asiatica*, which, according to the rickettsial classification guidelines (OrthoANI > 99.19%), indicate a clear species association toward *R. helvetica* (78).

From the eight positive specimens, we reconstructed a complete genome representation of *R. helvetica*. Additionally, two genomes (DK2 and OB144) were obtained using Vero cell lines. The 10 *de novo* assemblies exhibited high completeness, with BUSCO scores exceeding 98.9%. Genome sizes ranged from 1,372,202 bp to 1,372,575 bp, approximately 2.2–2.5 kb larger than the C9P9 reference genome (1,369,927 bp). Overall coverage ranged from 37,40× to 493,90× (File S1.10). All new *R. helvetica* genomes contained 3 rRNAs, 33 tRNAs, and between 1,630 and 1,664 protein-coding genes per assembly, including a complete set of genes for *de novo* folate biosynthesis (File S1.4). The overall coverage of these regions was analyzed, and no deviations were observed, indicating that an assembly error is unlikely.

Consistent with the reference strain, all assemblies included a plasmid of 47184–48,139 bp, encoding 56–59 proteins, including 39–40 hypothetical proteins. Plasmid association was verified using the reference strain and the DK2 and OB144 strains, all of which were obtained from culturing in Vero cell lines. Overall ANI values between the strains, including the plasmids, were on average 99.9% (File S1.12). No plasmid-specific reads were detected in the *R. helvetica*-negative specimens, indicating an association between *R. helvetica* and the plasmids identified.

Genome variant calling among the eight *l. ricinus*-derived genome assemblies revealed a total number of variants ranging from 63 to 82 (540 in total). Of these, 155 single nucleotide polymorphisms (SNPs), 344 insertions, 19 deletions, and 18 complex variants were detected (Table 2).

Overall, the new (meta)genome assemblies revealed a very high level of genomic collinearity, among themselves, and with the *R. helvetica* reference strain (Fig. 3 panel C). Most variation is due to what appears to be gene-splitting events, leading to two or more open reading frames. Between the eight genome assemblies that were directly obtained from the *l. ricinus* specimens, we found three such gene-splitting events. In contrast, we observed a notable increase of such events in the genomes of the DK2 strain and the C9P9 reference strain both obtained using infected Vero cells (File S1.6).

Due to such gene splitting events, the DK2 genome encodes a non-functional RickA protein (File S1.7). RickA protein plays a role in actin-based cell-to-cell motility, and *rickA* is therefore considered to be a putative virulence gene (34, 79). Further inspection of the new *R. helvetica* genomes revealed high genetic variability in the *rickA* gene due to copy number variation of a 33-nt-long repeat motif encoding 11 amino acids including five adjacent prolines ranging in number from 5 (Ir_dr4) to 13 (OB144).

To gain more insight into this genetic variability, we sequenced the fragment of the *rickA* gene containing the repeat motif from 61 questing nymphs gathered from
<table><caption>TABLE 2 Genetic variation observed in 11 genomes of *M. mitochondrii* and 8 genomes of *R. helvetica* analyzed in this study</caption><thead><tr><th>Strain</th><th>Complex<sup>a</sup></th><th>Deletions</th><th>Insertions</th><th>SNPs</th><th>Variant total<sup>b</sup></th><th>Genome size</th><th>% difference</th></tr></thead><tbody><tr><td><em>M. mitochondrii</em></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>F2</td><td>25</td><td>170</td><td>92</td><td>317</td><td>604</td><td>1,190,761</td><td>0.0059</td></tr><tr><td>Ir_d2</td><td>29</td><td>186</td><td>112</td><td>606</td><td>934</td><td>1,268,036</td><td>0.0664</td></tr><tr><td>Ir_d3</td><td>38</td><td>186</td><td>117</td><td>570</td><td>911</td><td>1,191,611</td><td>0.0066</td></tr><tr><td>Ir_d4</td><td>35</td><td>187</td><td>112</td><td>625</td><td>962</td><td>1,190,364</td><td>0.0055</td></tr><tr><td>Ir_d7</td><td>30</td><td>183</td><td>102</td><td>362</td><td>680</td><td>1,190,694</td><td>0.0058</td></tr><tr><td>Ir_d9</td><td>36</td><td>185</td><td>98</td><td>368</td><td>687</td><td>1,190,746</td><td>0.0058</td></tr><tr><td>Ir_f1</td><td>26</td><td>173</td><td>95</td><td>340</td><td>635</td><td>1,190,725</td><td>0.0058</td></tr><tr><td>Ir_f3</td><td>22</td><td>170</td><td>92</td><td>317</td><td>601</td><td>1,190,732</td><td>0.0058</td></tr><tr><td>Ir_f6</td><td>33</td><td>181</td><td>110</td><td>559</td><td>883</td><td>1,190,835</td><td>0.0059</td></tr><tr><td>Ir_f11</td><td>28</td><td>177</td><td>100</td><td>343</td><td>648</td><td>1,190,674</td><td>0.0058</td></tr><tr><td>Irf_16</td><td>33</td><td>191</td><td>116</td><td>622</td><td>964</td><td>1,190,312</td><td>0.0055</td></tr><tr><td><em>R. helvetica</em></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>Ir_d1</td><td>3</td><td>2</td><td>43</td><td>17</td><td>65</td><td>1,372,533</td><td>0.0018</td></tr><tr><td>Ir_d4</td><td>2</td><td>2</td><td>43</td><td>15</td><td>63</td><td>1,372,384</td><td>0.0017</td></tr><tr><td>Ir_d5</td><td>2</td><td>2</td><td>41</td><td>19</td><td>64</td><td>1,372,519</td><td>0.0018</td></tr><tr><td>Ir_d9</td><td>3</td><td>2</td><td>43</td><td>18</td><td>66</td><td>1,372,522</td><td>0.0018</td></tr><tr><td>Ir_f1</td><td>2</td><td>3</td><td>44</td><td>21</td><td>71</td><td>1,372,458</td><td>0.0018</td></tr><tr><td>Ir_f11</td><td>2</td><td>2</td><td>43</td><td>13</td><td>60</td><td>1,372,442</td><td>0.0018</td></tr><tr><td>Ir_f12</td><td>2</td><td>2</td><td>43</td><td>21</td><td>69</td><td>1,372,202</td><td>0.0016</td></tr><tr><td>Ir_f16</td><td>2</td><td>4</td><td>44</td><td>31</td><td>82</td><td>1,372,575</td><td>0.0019</td></tr></tbody></table>

various forest and dune regions (File S1). The results confirmed the observed copy number variation, which varied from 5 to 11. No significant compositional difference was observed between dune and forest areas (F 1,7 = 1.2, P = 0.3772), as well as between eight sampled locations with varying *R. helvetica* prevalence (F 1,7 = 0.4, P = 0.9101).

## Genome-scale metabolic modeling of *M. mitochondrii* and *R. helvetica* reveals metabolic dependencies

Due to their adaptation to rely on the host organism for survival, *M. mitochondrii* and *R. helvetica* have reduced genomes encoding reduced metabolic networks. The automatically generated GEM of *M. mitochondrii* (*Ir_d9_Mm*) is comprised of 1,126 reactions, 881 metabolites, and 266 genes, whereas the GEM of *R. helvetica* (*Ir_d9_Rh*) consists of 1,123 reactions, 851 metabolites, and 274 genes (Table 3).

MEMOTE analysis (66) revealed that both GEMs achieved high-quality overall scores of 84% (*Ir_d9_Mm*) and 68% (*Ir_d9_Rh*), respectively. (File S2). The network motifs of the two models point to differences in central carbon and fatty acid metabolisms. For instance, the central carbon metabolism appears to be more tightly connected to nitrogen metabolism in the *Ir_d9_Rh* GEM compared with *Ir_d9_Mm* GEM. Despite these differences, the two GEMs share 802 reactions. Each GEM also contained unique reactions and pathways: 300 for *Ir_d9_Rh* GEM and 429 for *Ir_d9_Mm* GEM. More specifically, the *Ir_d9_Mm* GEM contained unique reactions and pathways related to pantothenic acid (vitamin B5) and niacin (vitamin B3). Reaction essentiality analysis indicates that 8.70% and 8.17% of the reactions in *M. mitochondrii* and *R. helvetica*, respectively, are essential for *in silico* growth (File S2) and suggests that reactions linked to lysine synthesis pathways are vital for the growth of *R. helvetica* but not for *M. mitochondrii* (Fig. 4). Our analysis revealed that in addition to lysine, *R. helvetica* depends on the host for several amino acids such as valine, tyrosine, and isoleucine, reflecting its extensive auxotrophy (File S2).

<footnote><sup>a</sup>Complex: the number of multiple nucleotide polymorphisms.</footnote>

<footnote><sup>b</sup>Variant total: the total amount of variants between the newly obtained and reference genome.</footnote>
<table><caption>TABLE 3 Summary of genome-scale metabolic models of *M. mitochondrii* and *R. helvetica*</caption><thead><tr><th>Parameter</th><th>Ir_d9_Mm</th><th>Ir_d9_Rh</th></tr></thead><tbody><tr><td>Total no. of reactions</td><td>1126</td><td>1123</td></tr><tr><td>Total no. of metabolites</td><td>881</td><td>851</td></tr><tr><td>Total no. of genes</td><td>266</td><td>274</td></tr><tr><td>Reactions</td><td></td><td></td></tr><tr><td>Transport</td><td>349</td><td>399</td></tr><tr><td>Metabolic</td><td>630</td><td>581</td></tr><tr><td>Without GPR associations</td><td>337</td><td>353</td></tr><tr><td>Metabolites</td><td></td><td></td></tr><tr><td>Unique metabolites</td><td>881</td><td>851</td></tr><tr><td>Medium components</td><td>144</td><td>139</td></tr><tr><td>Genes: enzyme complexes</td><td>36</td><td>49</td></tr></tbody></table>

## Cophylogenetic analysis between *I. ricinus* and *symbionts*

The phylogenetic analysis of the mitogenomic data generated for 16 *I. ricinus* ticks was performed, including 33 mitochondrial genomes of *I. ricinus* obtained from the public databases. The results suggest a lack of a locality-dependent structure based on habitat type (dune and forest) from where ticks were collected (File S4; Fig. S4.1). The phylogenetic tree, however, did show an evolutionary distance between the different reconstructed mitochondrial genomes from the different samples.

To determine whether symbiont-tick pairs are co-evolved, the eMPress phylogenetic tree alignment software, following the duplication-transfer-loss model, was used

<figure><img src="image_5.png" /></figure>

<figure><img src="image_6.png" /><figcaption>FIG 4 Genome-scale metabolic constructions and analysis of *M. mitochondrii* and *R. helvetica*. Heatmap comparing *M. mitochondrii* (Ir_d9_Mm GEM) and *R. helvetica* (Ir9_Rh GEM) where a subset of the essential reactions is depicted. Furthermore, the blue and red panels represent which reactions are not essential or essential, respectively (A). Note: complete heatmap comparing the reaction essentiality of the two GEMs is contained in the supplementary material. Network map showing the network map of lysine synthesis pathway highlighting specific essential reactions in the Ir_d9_Mm GEM for *M. mitochondrii* (Ir_d9_Mm GEM) (B).</figcaption></figure>
to compare the phylogenetic relationships of symbionts and ticks. This analysis was preceded by the reconstruction of phylogenetic trees for *M. mitochondrii* and *R. helvetica* genomes obtained in this study (see File S4; Fig. S4.1 and S4.2). The reconciliation analysis was performed separately for two pairs: *M. mitochondrii*–*I. ricinus* and *R. helvetica*–*I. ricinus* (Fig. 5).

The results showed that a total of 12 events, including eight cospeciation events, two transfers, and two losses, occurred between *M. mitochondrii* and *I. ricinus* over evolutionary time; of these events, cospeciation was robustly supported in six out of eight events. The results of this analysis refuted the null hypothesis that the tick and *M. mitochondrii* trees and tip associations are formed due to chance at the 0.01 level in host-symbiont relationships (*P*-value = 0.0099). Therefore, we concluded that these *I. ricinus*–*M. mitochondrii* pairs have coevolved. The results for the *R. helvetica*–*I. ricinus* pair suggest rejection of the coevolution hypothesis (*P*-value = 0.77).

<figure><img src="image_7.png" /><figcaption>**FIG 5** Cophylogenetic analysis between *I. ricinus* ticks and its symbionts. (A) Phylogenetic reconciliation of *I. ricinus* (displayed in black) and *M. mitochondrii* (in blue). (B) Phylogenetic reconciliation of *I. ricinus* (displayed in black) and *R. helvetica* (in blue); (C and D) *P*-value histogram of host and *M. mitochondrii*; (E and F) *P*-value histogram of host and *R. helvetica*; the optimal reconciliation cost of the coevolution trees is indicated with a red line, and the optimal cost of the same trees constructed with tip associations permuted at random is shown in blue columns.</figcaption></figure>
## Rickettsia helvetica localization in adult female ticks

The distribution of *R. helvetica* in adult female ticks was studied to examine the nature and uniformity of its presence in tick tissues and to identify possible transmission routes. *R. helvetica* was successfully visualized using laser microscopy in five out of the 15 females (Fig. 6). These ticks were also qPCR positive for *R. helvetica*. All controls were negative (File S1.8). *R. helvetica* was observed in salivary glands, Malpighian tubules, the gut, and reproductive organs (File S5). The distribution of *R. helvetica* in tick organs varied between specimens. In the reproductive system, *R. helvetica* was observed in the female oviduct (two ticks, Fig. 6C and D) and in the *receptaculum seminis* (RS) (three ticks, Fig. 6A and B), which acts as the storage site for sperm introduced during copulation.

## DISCUSSION

### Endosymbiont genome sequencing

This study explores the transmission dynamics and lifestyle of symbionts associated with *I. ricinus*. Although traditional approaches for whole-genome sequencing of obligate intracellular bacteria have often relied on isolation and host-based propagation, recent advances have enabled direct sequencing without these steps (80, 81). In this study, we applied a high-depth, long-read in combination with a short-read direct deep-sequencing strategy to achieve the closed assembly of a relatively large set of samples. This approach allows for a more efficient reconstruction of the *M. mitochondrii* and *R. helvetica* genomes, providing insights into their transmission dynamics and lifestyle.

In this study, *M. mitochondrii* was detected in 15 specimens, except for a single specimen molted under laboratory conditions. This absence aligns with previous observations that both the prevalence and abundance of *M. mitochondrii* decrease in laboratory-maintained tick colonies (82–84). Such reductions suggest that the bacterium plays a role specific to natural environments, potentially linked to host interactions, environmental pressures, or nutrient availability absent in lab settings.

Previous research indicates that laboratory interventions can significantly diminish *M. mitochondrii* populations. For instance, a study by Guizzo et al. demonstrated that administering tetracycline to *I. ricinus* ticks via an *ex vivo* membrane feeding system significantly reduced the abundance of *M. mitochondrii* in their ovaries (82). This reduction persisted across subsequent generations, leading to a complete elimination of the symbiont. Similarly, Ninio et al. observed that antibiotic treatment of *I. ricinus* resulted in a decreased load of *M. mitochondrii* following a blood meal (85). These findings suggest that laboratory interventions, such as antibiotic administration, can effectively reduce or eliminate *M. mitochondrii* in tick populations.

In addition to *M. mitochondrii*, *R. helvetica* was present in eight specimens, whereas *Rickettsiella* species were prominently present in three specimens. A high overall abundance of *Pseudomonas*, *Mycobacterium*, and *Mycobacteroides* species was also observed. It is important to note that the sensitivity of a Kraken2 analysis is significantly influenced by the type of genetic data available in the nt database—whether 16S sequences, individual gene fragments, or complete genomes—indicating that Kraken read counts cannot be directly used to estimate relative species abundances.

As there are reference genomes available for *M. mitochondrii* and *R. helvetica*, employing fixed-length Illumina read pairs as input allowed us to directly utilize the Kraken2 output for estimating whether there exists adequate short-read genome coverage to reconstruct the complete microbial genomes of interest. In alignment with the outcomes from Kraken2, we were able to obtain 11 complete and four draft genomes of *M. mitochondrii* and, coinciding with qPCR analysis results, eight complete *R. helvetica* genomes along with their associated plasmids.

### M. mitochondrii

The genomes of *M. mitochondrii* as presented in this study are approximately 5% larger and have a structural arrangement that differs from the published assembly (accession
<figure><img src="image_8.png" /><figcaption>**FIG 6** Localization of *R. helvetica* in adult females. Blue: DAPI, green dots: Rickettsia-specific binding. (A) *R. helvetica*-positive *receptaculum seminis* X40 A' *R. helvetica*-positive *receptaculum seminis* ×100, (B) control *receptaculum seminis*, (C) *R. helvetica*-positive oviduct X40, and (D) control oviduct. Note that tissue-specific autofluorescence can be observed in the treatment and control groups. Scale bar: 10 µm</figcaption></figure>

no.: NC_015722). The observed differences in genome size and structure could be due to the use of different technological platforms, assembly methods, and geographic origin (the Netherlands vs. Italy). The extraction of 11 complete *M. mitochondrii* genomes, each from a single tick, provides insight into the genomic potential and structural variability of this symbiotic bacterium in its tick host.

*M. mitochondrii* exhibits a unique feature to colonize multiple compartments in a tick body, including ovaries, salivary glands (86), Malpighian tubules, tracheae, and guts (13). In this context, each genome generated in our study represents a mixture of genomes from bacteria inhabiting different ecological niches within the tick body. Subsequently, in the eMPRess reconciliation analysis, we showed a significant coevolutionary relationship between *M. mitochondrii* and its tick host. This finding supports the notion of a long-standing symbiotic relationship between *M. mitochondrii* and *I. ricinus*, suggesting that vertical transmission is likely the primary mode of symbiont transfer.
The strong cophylogeny observed between *M. mitochondrii* and its host's mitochondria in this study further supports the idea that vertical transmission is the dominant mode of inheritance, opposing the hypothesis of horizontal transmission.

ANI analysis of the *M. mitochondrii* complete draft genomes revealed two main groups of genome arrangements: one consisting of *M. mitochondrii* isolated from *Ir_f1*, *Ir_d3*, *Ir_d7*, *Ir_f3*, and *Ir_f11*, and another consisting of *Ir_d2*, *Ir_d4*, and *Ir_f16* (File S1.12). This grouping further supported the strong cophylogeny observed between *M. mitochondrii* and its host's mitochondria (File S4.1).

Although *M. mitochondrii* has been detected in tracheae, salivary glands, and rostrum, as well as in the blood of experimentally infested rabbits, horizontal transmission is unlikely to play a major role in its spread within natural tick populations. However, given the localized scope of this data set, broader conclusions regarding the absence of horizontal transmission would require further investigation across diverse populations and environments.

## R. helvetica

The prevalence of *R. helvetica* in questing ticks has been shown to greatly differ between dune and forest areas in the Netherlands (27, 33, 87). Although our findings point to some heterogeneity in the *rickA* virulence gene, there seems to be no genetic basis for the differences in the prevalence of *R. helvetica* in dune versus forest areas. Thus, two remaining possibilities should be further explored. First, the possibility of *R. helvetica* giving ticks in dune areas an adaptive advantage in relation to local abiotic conditions has not been fully explored. The role of *R. helvetica* in this case could be related to nutrient acquisition, waste elimination, or physical characteristics. For example, *Rickettsia* has been shown to enhance tick motility in *Dermacentor variabilis* (88). Second, the possibility of interactions with other symbionts within the tick host should be explored as *R. helvetica* has been shown to modulate the microbiome of *I. ricinus* (89), and it has been found to be negatively correlated with certain symbionts such as *Spiroplasma ixodetis* and *M. mitochondrii* (2, 33). It has been hypothesized that perhaps the biological roles of *S. ixodetis* and *R. helvetica* overlap, making their simultaneous maintenance energetically superfluous for the tick host (33), but this theory has yet to be proven. Nor has the option of other bacterial interactions been fully researched.

Our comparative analyses showed near absolute genetic homogeneity among the generated *R. helvetica* genomes (99.96%–100%), and only a slight deviation from the published strain C9P9. In the latter case, the observed differences in genome size could be due to different isolation and assembly methods or even geographic origin. The congruous genetic build-up among the *R. helvetica* genomes was punctuated almost exclusively by the gene encoding for the surface protein RickA. The sequence of the *rickA* gene showed high genetic variation between the genomes generated in this study as well as between genes obtained from 61 additional ticks from the Netherlands. The RickA protein has been shown to induce actin polymerization and hence actin-based motility in some *Rickettsia* species (37, 80). However, a previous study found the RickA protein in *R. helvetica* isolate AS819 to be insufficient to promote actin-based motility in vitro, which spread instead through cell breakdown and bacterial release (26). The inability of the studied isolate to spread through actin-based motility was explained by a disrupted *sca2* gene, another virulence gene believed to have a more prominent role than the *rickA* gene in the cell-to-cell spread of *Rickettsias* (90). Our genomes revealed no truncated *sca2* sequences; however, the *rickA* gene greatly differed in length due to a repeat motif. We hypothesize these repeats may affect the functionality of the RickA protein within *R. helvetica*, giving rise to variation in virulence. This observation is supported by the dose-dependent manner in which RickA induced cell motility in *R. conorii* (37). Moreover, when expressing a heterologous *rickA* gene, *R. bellii* moved at higher velocities than its wild-type counterparts (91). If the diverse *rickA* sequence length translates into different motility-associated phenotypes, it could explain the seemingly sporadic pathogenicity we observe for *R. helvetica*, as well as its low and inconsistent prevalence rate in wildlife
hosts. This heterogeneity could also explain why two different mouse models showed no substantial *rickettsia* following an infection challenge with *R. helvetica* (92) and should be considered when exploring the infection potential of *R. helvetica* strains.

To date, no wildlife reservoir has been found for *R. helvetica*. This has led to the assumption that *I. ricinus*, which transmits *R. helvetica* transovarially (93), acts as its main reservoir in nature (27). This close tick-symbiont interaction was not reflected in our findings, which showed the absence of co-phylogeny between *R. helvetica* and the mitochondrial genomes of *I. ricinus*. The lack of co-evolution between *I. ricinus* and *R. helvetica* may be partly due to a mixed transmission strategy encompassing vertical and horizontal routes. The presence of *R. helvetica* in the RS, as shown by our immunofluorescence findings, strongly suggests the possibility of paternal transmission for *R. helvetica* in *I. ricinus*. Paternal transmission of *Rickettsia* has been identified in the leafhopper *Nephotettix cincticeps* (94), and it has been hypothesized that a mixed mode of transmission (maternal and parental) might allow for the mixing of different symbiont lineages, thus allowing for the evolution of virulent variants (95). *Rickettsia* has also been detected in the spermatids of *Dermacentor silvarum* and *I. ricinus* ticks (96, 97). In this study, adult females were studied, as pathogens are thought to be transmitted to humans mainly by nymphs and female ticks. The possibility of a different genetic makeup for *Rickettsia* populations found in reproductive tick organs should be considered in future studies in both male and female ticks.

## Symbiont metabolic integration

Genome-scale metabolic models (GEMs) provide a means to uncover new biological insights that go beyond what genomics can offer. A key application involves predicting the viability of an organism under specific conditions. This simulation approach has been utilized to identify potential metabolic drug targets, enabling the effective targeting of pathogens (98).

Due to their long-term relationship with the tick host, both *M. mitochondrii* and *R. helvetica* underwent genome reduction, leading to incomplete metabolic networks. The analysis suggests that *M. mitochondrii's* *in silico* growth necessitates the inclusion of pantothenate (vitamin B5) and glycerol-3-phosphate (G3P) in the growth medium. G3P, a key metabolite in metabolic pathways such as glycolysis, gluconeogenesis, and lipid metabolism, is abundantly available in mitochondria. Its role is pivotal in the glycerol phosphate shuttle, facilitating the transfer of reducing equivalents to the mitochondrial matrix. Pantothenate is vital for energy metabolism and the synthesis of coenzyme A (CoA), a coenzyme integral to mitochondrial metabolic reactions. Our findings about the growth requirements of *M. mitochondrii* are in line with the evolutionary adaptation of this symbiont to reside in its host's mitochondria. Nevertheless, as ticks cannot synthesize pantothenate, there could be potentially other bacterial symbionts producing this vitamin to be utilized by *M. mitochondrii*. The GEM for *R. helvetica* highlights a dependency on lysine for *in silico* growth, as demonstrated in Fig. 4B. This requirement likely results from an incomplete lysine biosynthesis pathway, a characteristic shared with *Rickettsia prowazekii* (99, 100), as previously identified. Research has consistently shown that *Rickettsia* species lack comprehensive amino acid biosynthesis pathways, particularly for essential amino acids such as lysine (72, 101). Therefore, *R. helvetica* must obtain lysine from the blood meal ingested by its tick host or from other tick symbionts; however, the latter has not been explored. In some cases, *Rickettsia* species have been found to carry out only the initial stages of lysine biosynthesis, leading to the production of diaminopimelate rather than lysine itself, necessitating the acquisition of lysine from their external environment (101).

Simultaneously, we examined how these symbionts could contribute to their tick host by providing essential micronutrients. All *R. helvetica* genomes harbor a complete set of genes for *de novo* folate biosynthesis, in line with previous findings (13, 102). *M. mitochondrii* encodes a (near) complete folate biosynthetic pathway and a likely complete biotin biosynthetic pathway.
## Conclusion

Ticks harbor and transmit specialized symbiotic microorganisms that are crucial for their hematophagous lifestyle and occasionally serve as causal agents of diseases in humans and animals. Although these symbionts are challenging to culture and manipulate ex vivo, our study shows how deep sequencing allows for the direct extraction of complete genomes from tick specimens collected from their natural habitats, providing fresh insights into their structural genome variation, host interactions, and modes of transmission.

## ACKNOWLEDGMENTS

B.N., W.T.S., P.J.S., and J.J.K. acknowledge the Dutch Research Council (NWO) and Wageningen University & Research for their financial contribution to the Unlock initiative (NWO: 184.035.007). T.A. and H.S. were supported by funding from the Dutch Ministry of Health, Welfare, and Sport.

## AUTHOR AFFILIATIONS

¹Department of Veterinary Sciences, Faculty of Agrobiology, Food and Natural Resources, Czech University of Life Sciences, Prague, Prague, Czechia

²NVWA-Centrum Monitoring Vectoren, Wageningen, the Netherlands

³Centre for Infectious Diseases (Clb), National Institute for Public Health and the Environment (RIVM), Bilthoven, Utrecht, the Netherlands

⁴Laboratory of Systems and Synthetic Biology, Wageningen University & Research, Wageningen, Gelderland, the Netherlands

⁵UNLOCK, Wageningen University & Research, Wageningen, Gelderland, the Netherlands

⁶Future Genomics Technologies BV, Leiden, the Netherlands

⁷Laboratoire de Santé Animale, École Nationale Vétérinaire d'Alfort, UMR BIPAR, Anses INRAE, Maisons-Alfort, France

⁸Bundeswehr Institute of Microbiology, Munich, Bavaria, Germany

## AUTHOR ORCIDs

Hein Sprong http://orcid.org/0000-0002-0218-4320

Jasper J. Koehorst http://orcid.org/0000-0001-8172-8981

## FUNDING

<table><thead><tr><th>Funder</th><th>Grant(s)</th><th>Author(s)</th></tr></thead><tbody><tr><td>Nederlandse Organisatie voor Wetenschappelijk Onderzoek</td><td>184.035.007</td><td>William T. Scott Jr.<br>Bart Nijsse<br>Peter J. Schaap<br>Jasper J. Koehorst</td></tr><tr><td>Ministerie van Volksgezondheid, Welzijn en Sport</td><td></td><td>Tal Azagi<br>Hein Sprong</td></tr></tbody></table>

## AUTHOR CONTRIBUTIONS

Paulina M. Lesiczka, Conceptualization, Data curation, Formal analysis, Investigation, Methodology, Validation, Visualization, Writing – original draft, Writing – review and editing | Tal Azagi, Conceptualization, Data curation, Formal analysis, Investigation, Methodology, Supervision, Writing – original draft, Writing – review and editing | Aleksandra I. Krawczyk, Conceptualization, Methodology, Project administration, Writing – review and editing | William T. Scott, Jr., Investigation, Methodology, Resources, Software, Validation, Writing – original draft | Ron P. Dirks, Data curation, Formal analysis, Resources | Gerhard Dobler, Conceptualization, Funding acquisition, Investigation,
Methodology, Project administration | Bart Nijsse, Data curation, Formal analysis, Investigation, Methodology, Resources, Software, Validation, Visualization, Writing – original draft, Writing – review and editing | Hein Sprong, Conceptualization, Funding acquisition, Investigation, Project administration, Supervision, Writing – original draft, Writing – review and editing.

## DATA AVAILABILITY

The data supporting the findings of this study have been deposited in the European Nucleotide Archive (ENA). This data set includes Oxford Nanopore Technologies (ONT) and Illumina sequencing data for multiple *Ixodes ricinus* specimens. The project number is PRJEB70514. ONT accession numbers are as follows: ERR12319604, ERR12319605, ERR12319607, ERR12319608, ERR12319610, ERR12319612, ERR12319614, ERR12319616, ERR12319618, ERR12319620, ERR12319623, ERR12319624, ERR12319627, ERR12319628, ERR12319631, ERR12319632, ERR12319593, ERR12319595, ERR12319597, ERR12319598, ERR12319599, ERR12319600, ERR12319602. Illumina accession numbers are as follows: ERR12319603, ERR12319606. ERR12319609, ERR12319611, ERR12319615, ERR12319617, ERR12319619, ERR12319621, ERR12319625, ERR12319626, ERR12319629, ERR12319630, ERR12319633, ERR12319622, ERR12319594, ERR12319596, ERR12319601. *Midichloria mitochondrii* accession numbers are as follows: GCA_963853845, GCA_963854075, GCA_963853965, GCA_963853945, GCA_963853935, GCA_963854045, GCA_963854055, GCA_963853925, GCA_963853985, GCA_963854105. *Rickettsia helvetica* accession numbers are as follows: GCA_963854155, GCA_963853705, GCA_963853725, GCA_963853685, GCA_963853695, GCA_963853745, GCA_963853735, GCA_963853715, GCA_963854155, GCA_963970025. The code used for the model simulations can be found at https://git.wur.nl/unlock/projects/ticks, with supplemental text and figures at https://git.wur.nl/unlock/projects/ticks/-/tree/main/supplementary_material. The GEMs as well as accompanying MEMOTE (56) and FROG analysis (11) reports can be found in a public repository on BioModels: https://www.ebi.ac.uk/biomodels/. For the *M. mitochondrii* GEM, see here: MODEL2404170001. For the *R. helvetica* GEM, see MODEL2404170002.

## ADDITIONAL FILES

The following material is available online.

### Supplemental Material

File S1 (mSystems00507-25-s0001.xlsx). Information about the strains; coverage; BUSCO, ANI, and TERA analysis; genome accessions; eggNOG annotation; RickA; and co-linearity analysis.

File S2 (mSystems00507-25-s0002.pdf). Genome-scale modeling details.

File S3 (mSystems00507-25-s0003.pdf). fastANI visualization.

File S4 (mSystems00507-25-s0004.pdf). Phylogenetic tree of the mitochondrial *Ixodes ricinus*.

File S5 (mSystems00507-25-s0005.pdf). *Rickettsia helvetica* localization in adult female ticks.

## REFERENCES

1. Sprong H, Azagi T, Hoornstra D, Nijhof AM, Knorr S, Baarsma ME, Hovius JW. 2018. Control of Lyme borreliosis and other *Ixodes ricinus*-borne diseases. Parasit Vectors 11:145. https://doi.org/10.1186/s13071-018-2744-5

2. Krawczyk AI, Röttjers S, Coimbra-Dores MJ, Heylen D, Fonville M, Takken W, Faust K, Sprong H. 2022. Tick microbial associations at the crossroad of horizontal and vertical transmission pathways. Parasit Vectors 15:380. https://doi.org/10.1186/s13071-022-05519-w

3. Azagi T, Hoornstra D, Kremer K, Hovius JWR, Sprong H. 2020. Evaluation of disease causality of rare *Ixodes ricinus*-borne infections in Europe. Pathogens 9:E150. https://doi.org/10.3390/pathogens9020150

4. Tijsse-Klasen E, Koopmans MPG, Sprong H. 2014. Tick-borne pathogen – reversed and conventional discovery of disease. Front Public Health 2:73. https://doi.org/10.3389/fpubh.2014.00073

5. Benelli G. 2020. Pathogens manipulating tick behavior-through a glass, darkly. Pathogens 9:664. https://doi.org/10.3390/pathogens9080664

6. Bonnet SI, Binetruy F, Hernández-Jarguin AM, Duron O. 2017. The tick microbiome: why non-pathogenic microorganisms matter in tick
biology and pathogen transmission. Front Cell Infect Microbiol 7:236. https://doi.org/10.3389/fcimb.2017.00236

7. Duron O. 2024. Nutritional symbiosis in ticks: singularities of the genus *Ixodes*. Trends Parasitol 40:696-706. https://doi.org/10.1016/j.pt.2024.06.006

8. Lewis D. 1979. The detection of rickettsia-like microorganisms within the ovaries of female *Ixodes ricinus* ticks. Z Parasitenkd 59:295-298. http://doi.org/10.1007/BF00927523

9. Sacchi L, Bigliardi E, Corona S, Beninati T, Lo N, Franceschi A. 2004. A symbiont of the tick *Ixodes ricinus* invades and consumes mitochondria in a mode similar to that of the parasitic bacterium *Bdellovibrio bacteriovorus*. Tissue Cell 36:43-53. https://doi.org/10.1016/j.tice.2003.08.004

10. Sassera D, Beninati T, Bandi C, Bouman EAP, Sacchi L, Fabbi M, Lo N. 2006. "Candidatus Midichloria mitochondrii", an endosymbiont of the tick *Ixodes ricinus* with a unique intramitochondrial lifestyle. Int J Syst Evol Microbiol 56:2535-2540. https://doi.org/10.1099/ijs.0.64386-0

11. Sassera D, Lo N, Bouman EAP, Epis S, Mortarino M, Bandi C. 2008. "Candidatus Midichloria" endosymbionts bloom after the blood meal of the host, the hard tick *Ixodes ricinus*. Appl Environ Microbiol 74:6138-6140. https://doi.org/10.1128/AEM.00248-08

12. Floriano AM, Batisti Bifignandi G, Castelli M, Olivieri E, Clementi E, Comandatore F, Rinaldi L, Opara M, Plantard O, Palomar AM, Noël V, Vijay A, Lo N, Makepeace BL, Duron O, Jex A, Guy L, Sassera D. 2023. The evolution of intramitochondriality in *Midichloria* bacteria. Environ Microbiol 25:2102-2117. https://doi.org/10.1111/1462-2920.16446

13. Olivieri E, Epis S, Castelli M, Varotto Boccazzi I, Romeo C, Desirò A, Bazzocchi C, Bandi C, Sassera D. 2019. Tissue tropism and metabolic pathways of *Midichloria mitochondrii* suggest tissue-specific functions in the symbiosis with *Ixodes ricinus*. Ticks Tick Borne Dis 10:1070-1077. https://doi.org/10.1016/j.ttbdis.2019.05.019

14. Bazzocchi C, Mariconti M, Sassera D, Rinaldi L, Martin E, Cringoli G, Urbanelli S, Genchi C, Bandi C, Epis S. 2013. Molecular and serological evidence for the circulation of the tick symbiont *Midichloria* (Rickettsiales: Midichloriaceae) in different mammalian species. Parasit Vectors 6:350. https://doi.org/10.1186/1756-3305-6-350

15. Murray GGR, Weinert LA, Rhule EL, Welch JJ. 2016. The phylogeny of *Rickettsia* using different evolutionary signatures: how tree-like is bacterial evolution? Syst Biol 65:265-279. https://doi.org/10.1093/sysbi o/syv084

16. Nilsson K., Lindquist O, Pålsson C. 1999. Association of *Rickettsia helvetica* with chronic perimyocarditis in sudden cardiac death. The Lancet 354:1169-1173. https://doi.org/10.1016/S0140-6736(99)04093-3

17. Fournier PE, Grunnenberger F, Jaulhac B, Gastinger G, Raoult D. 2000. Evidence of *Rickettsia helvetica* infection in humans, eastern France. Emerg Infect Dis 6:389-392. https://doi.org/10.3201/eid0604.000412

18. Fournier P-E, Allombert C, Supputamongkol Y, Caruso G, Brouqui P, Raoult D. 2004. Aneruptive fever associated with antibodies to *Rickettsia helvetica* in Europe and Thailand. J Clin Microbiol 42:816-818. https://doi.org/10.1128/JCM.42.2.816-818.2004

19. Nilsson K. 2009. Septicaemia with *Rickettsia helvetica* in a patient with acute febrile illness, rash and myasthenia. J Infect 58:79-82. https://doi. org/10.1016/j.jinf.2008.06.005

20. Nilsson K, Elfving K, Pålsson C. 2010. *Rickettsia helvetica* in patient with meningitis, Sweden, 2006. Emerg Infect Dis 16:490-492. https://doi.org/ 10.3201/eid1603.090184

21. Nilsson K, Wallménius K, Pålsson C. 2011. Coinfection with *Rickettsia helvetica* and herpes simplex virus 2 in a young woman with meningoencephalitis. Case Rep Infect Dis 2011:469194. https://doi.org/10.115 5/2011/469194

22. Koetsveld J, Tijsse-Klasen E, Herremans T, Hovius JWR, Sprong H. 2016. Serological and molecular evidence for spotted fever group *Rickettsia* and *Borrelia burgdorferi sensu lato* co-infections in The Netherlands. Ticks Tick Borne Dis 7:371-377. https://doi.org/10.1016/j.ttbdis.2015.12. 010

23. Ocias LF, Dessau RB, Lebech A-M, Jørgensen CS, Petersen RF, Krogfelt KA. 2018. Evidence of rickettsiae in Danish patients tested for Lyme neuroborreliosis: a retrospective study of archival samples. BMC Infect Dis 18:325. https://doi.org/10.1186/s12879-018-3210-x

24. Esteves T, Aparicio G, Garcia-Patos V. 2016. Is there any association between sarcoidosis and infectious agents?: a systematic review and meta-analysis. BMC Pulm Med 16:165. https://doi.org/10.1186/s12890- 016-0332-z

25. Rizzoli A, Silaghi C, Obiegala A, Rudolf I, Hubálek Z, Földvári G, Plantard O, Vayssier-Taussat M, Bonnet S, Spitalská E, Kazimírová M. 2014. *Ixodes ricinus* and its transmitted pathogens in urban and peri-urban areas in Europe: new hazards and relevance for public health. Front Public Health 2:251. https://doi.org/10.3389/fpubh.2014.00251

26. Speck S, Kern T, Aistleitner K, Dilcher M, Dobler G, Essbauer S. 2018. *In vitro* studies of *Rickettsia*-host cell interactions: confocal laser scanning microscopy of *Rickettsia helvetica*-infected eukaryotic cell lines. PLoS Negl Trop Dis 12:e0006151. https://doi.org/10.1371/journal.pntd.00061 51

27. Sprong H, Wielinga PR, Fonville M, Reusken C, Brandenburg AH, Borgsteede F, Gaasenbeek C, van der Giessen JW. 2009. *Ixodes ricinus* ticks are reservoir hosts for *Rickettsia helvetica* and potentially carry fleas- borne *Rickettsia* species. Parasit Vectors 2:41. https://doi.org/10.1186/17 56-3305-2-41

28. Inokuma H, Seino N, Suzuki M, Kaji K, Takahashi H, Igota H, Inoue S. 2008. Detection of *Rickettsia helvetica* DNA from peripheral blood of Sika deer (*Cervus nippon yesoensis*) in Japan. J Wildl Dis 44:164-167. https://doi.org/10.7589/0090-3558-44.1.164

29. Boretti FS, Perreten A, Meli ML, Cattori V, Willi B, Wengi N, Hornok S, Honegger H, Hegglin D, Woelfel R, Reusch CE, Lutz H, Hofmann-Lehmann R. 2009. Molecular investigations of *Rickettsia helvetica* infection in dogs, foxes, humans, and *Ixodes* ticks. Appl Environ Microbiol 75:3230-3237. https://doi.org/10.1128/AEM.00220-09

30. Elfving K, Malmsten J, Dalin A-M, Nilsson K. 2015. Serologic and molecular prevalence of *Rickettsia helvetica* and *Anaplasma phagocytophilum* in wild cervids and domestic mammals in the central parts of Sweden. Vector Borne Zoonotic Dis 15:529-534. https://doi.org/10.108 9/vbz.2015.1768

31. Hornok S, Kováts D, Csörgő T, Meli ML, Gönczi E, Hadnagy Z, Takács N, Farkas R, Hofmann-Lehmann R. 2014. Birds as potential reservoirs of tick-borne pathogens: first evidence of bacteraemia with *Rickettsia helvetica*. Parasit Vectors 7:128. https://doi.org/10.1186/1756-3305-7-12 8

32. Michelet L, Delannoy S, Devillers E, Umhang G, Aspan A, Juremalm M, Chirico J, van der Wal FJ, Sprong H, Boye Pihl TP, Klitgaard K, Bødker R, Fach P, Moutailler S. 2014. High-throughput screening of tick-borne pathogens in Europe. Front Cell Infect Microbiol 4:103. https://doi.org/1 0.3389/fcimb.2014.00103

33. Krawczyk Al, Röttjers L, Fonville M, Takumi K, Takken W, Faust K, Sprong H. 2022. Quantitative microbial population study reveals geographical differences in bacterial symbionts of *Ixodes ricinus*. Microbiome 10:120. https://doi.org/10.1186/s40168-022-01276-1

34. Harris EK, Jirakanwisal K, Verhoeve VI, Fongsaran C, Suwanbongkot C, Welch MD, Macaluso KR. 2018. Role of Sca2 and RickA in the dissemination of *Rickettsia parkeri* in *Amblyomma maculatum*. Infect Immun 86:e00123-18. https://doi.org/10.1128/IAI.00123-18

35. Ireton K. 2013. Molecular mechanisms of cell-cell spread of intracellular bacterial pathogens. Open Biol 3:130079. https://doi.org/10.1098/rsob. 130079

36. Merhej V, Raoult D. 2011. Rickettsial evolution in the light of comparative genomics. Biol Rev Camb Philos Soc 86:379-405. https://doi.org/10. 1111/j.1469-185X.2010.00151.x

37. Gouin E, Egile C, Dehoux P, Villiers V, Adams J, Gertler F, Li R, Cossart P. 2004. The RickA protein of *Rickettsia conorii* activates the Arp2/3 complex. Nature 427:457-461. https://doi.org/10.1038/nature02318

38. Wang T, Zhang S, Pei T, Yu Z, Liu J. 2019. Tick mitochondrial genomes: structural characteristics and phylogenetic implications. Parasit Vectors 12:451. https://doi.org/10.1186/s13071-019-3705-3

39. Kneubehl AR, Muñoz-Leal S, Filatov S, de Klerk DG, Pienaar R, Lohmeyer KH, Bermúdez SE, Suriyamongkol T, Mali I, Kanduma E, Latif AA, Sarih M, Bouattour A, de León AAP, Teel PD, Labruna MB, Mans BJ, Lopez JE. 2022. Amplification and sequencing of entire tick mitochondrial genomes for a phylogenomic analysis. Sci Rep 12:19310. https://doi.org/10.1038/s41598-022-23393-5

40. Ahantarig A, Trinachartvanit W, Baimai V, Grubhoffer L. 2013. Hard ticks and their bacterial endosymbionts (or would be pathogens). Folia Microbiol (Praha) 58:419-428. https://doi.org/10.1007/s12223-013-022 2-1

41. Tian J, Hou X, Ge M, Xu H, Yu B, Liu J, Shao R, Holmes EC, Lei C, Shi M. 2022. The diversity and evolutionary relationships of ticks and tick-borne bacteria collected in China. Parasit Vectors 15:352. https://doi.org/10.1186/s13071-022-05485-3
42. Stavru F, Riemer J, Jex A, Sassera D. 2020. When bacteria meet mitochondria: The strange case of the tick symbiont *Midichloria mitochondrii*. Cell Microbiol 22:e13189. https://doi.org/10.1111/cmi.13189

43. Ankrah NYD, Luan J, Douglas AE. 2017. Cooperative metabolism in a three-partner insect-bacterial symbiosis revealed by metabolic modeling. J Bacteriol 199:e00872-16. https://doi.org/10.1128/JB.00872-16

44. Hillyard PD. 1996. Ticks of north-west Europe. Field Studies Council, Shrewsbury, England.

45. Krawczyk AI, Bakker JW, Koenraadt CJM, Fonville M, Takumi K, Sprong H, Demir S. 2020. Tripartite interactions among *Ixodaphagus hookeri*, *Ixodes ricinus* and deer: differential interference with transmission cycles of tick-borne pathogens. Pathogens 9:339. https://doi.org/10.3390/pathogens9050339

46. Lu J, Breitwieser FP, Thielen P, Salzberg SL. 2016. Bracken: estimating species abundance in metagenomics data. bioRxiv. https://doi.org/10.1101/051813

47. Milhano N, Palma M, Marcili A, Núncio MS, de Carvalho IL, de Sousa R. 2014. *Rickettsia lusitaniae* sp. nov. isolated from the soft tick *Ornithodoros erraticus* (Acarina: Argasidae). Comp Immunol Microbiol Infect Dis 37:189–193. https://doi.org/10.1016/j.cimid.2014.01.006

48. Ammerman NC, Beier-Sexton M, Azad AF. 2008. Laboratory maintenance of *Rickettsia rickettsii*. Curr Protoc Microbiol 11:3A.5.1–3A.5.21. https://doi.org/10.1002/9780471729259.mc03a05s11

49. J.J.Koehorst BN. 2024. (hybrid) metagenomics workflow. Available from: https:// workflowhub.eu/workflows/367

50. Goble C, Soiland-Reyes S, Bacall F, Owen S, Williams A, Eguinoa I, Droesbeke B, Leo S, Pireddu L, Rodriguez-Navas L, Fernández JM, Capella-Gutierrez S, Ménager H, Grüning B, Serrano-Solano B, Ewels P, Coppens F. 2021. Implementing FAIR digital objects in the EOSC-life workflow collaboratory. Zenodo. https://doi.org/10.5281/zenodo.4605654.

51. Andrews S. FastQC. Babraham Institute.

52. Chen S, Zhou Y, Chen Y, Gu J. 2018. fastp: an ultra-fast all-in-one FASTQ preprocessor. Bioinformatics 34:i884–i890. https://doi.org/10.1093/bioinformatics/bty560

53. Bushnell B. 2014. BBMap: a fast, accurate, splice-aware aligner

54. Li H. 2018. Minimap2: pairwise alignment for nucleotide sequences. Bioinformatics 34:3094–3100. https://doi.org/10.1093/bioinformatics/bty191

55. Hartemink N, van Vliet A, Sprong H, Jacobs F, Garcia-Martí I, Zurita-Milla R, Takken W. 2019. Temporal-spatial variation in questing tick activity in the Netherlands: the effect of climatic and habitat factors. Vector Borne Zoonotic Dis 19:494–505. https://doi.org/10.1089/vbz.2018.2369

56. Hansford KM, Fonville M, Gillingham EL, Coipan EC, Pietzsch ME, Krawczyk AI, Vaux AGC, Cull B, Sprong H, Medlock JM. 2017. Ticks and *Borrelia* in urban and peri-urban green space habitats in a city in southern England. Ticks Tick Borne Dis 8:353–361. https://doi.org/10.1016/j.ttbdis.2016.12.009

57. Kearse M, Moir R, Wilson A, Stones-Havas S, Cheung M, Sturrock S, Buxton S, Cooper A, Markowitz S, Duran C, Thierer T, Ashton B, Meintjes P, Drummond A. 2012. Geneious basic: an integrated and extendable desktop software platform for the organization and analysis of sequence data. Bioinformatics 28:1647–1649. https://doi.org/10.1093/bioinformatics/bts199

58. Thompson JD, Gibson TJ, Higgins DG. 2002. Multiple sequence alignment using ClustalW and ClustalX. Curr Protoc Bioinform 00:2.3.1–2.3.22 https://doi.org/10.1002/0471250953.bi0203s00

59. Minh BQ, Schmidt HA, Chernomor O, Schrempf D, Woodhams MD, von Haeseler A, Lanfear R. 2020. IQ-TREE 2: new models and efficient methods for phylogenetic inference in the genomic era. Mol Biol Evol 37:1530–1534. https://doi.org/10.1093/molbev/msaa015

60. Kalyaanamoorthy S, Minh BQ, Wong TKF, von Haeseler A, Jermiin LS. 2017. ModelFinder: fast model selection for accurate phylogenetic estimates. Nat Methods 14:587–589. https://doi.org/10.1038/nmeth.4285

61. Hoang DT, Chernomor O, von Haeseler A, Minh BQ, Vinh LS. 2018. UFBoot2: improving the ultrafast bootstrap approximation. Mol Biol Evol 35:518–522. https://doi.org/10.1093/molbev/msx281

62. Guindon S, Dufayard J-F, Lefort V, Anisimova M, Hordijk W, Gascuel O. 2010. New algorithms and methods to estimate maximum-likelihood phylogenies: assessing the performance of PhyML 3.0. Syst Biol 59:307–321. https://doi.org/10.1093/sysbio/syq010

63. Bah T. 2011. Inkscape: guide to a vector drawing program. Prentice Hall Press.

64. Santichaivekin S, Yang Q, Liu J, Mawhorter R, Jiang J, Wesley T, Wu Y-C, Libeskind-Hadas R. 2021. eMPress: a systematic cophylogeny reconciliation tool. Bioinformatics 37:2481–2482. https://doi.org/10.1093/bioinformatics/btaa978

65. Machado D, Andrejev S, Tramontano M, Patil KR. 2018. Fast automated reconstruction of genome-scale metabolic models for microbial species and communities. Nucleic Acids Res 46:7542–7553. https://doi.org/10.1093/nar/gky537

66. Lieven C, Beber ME, Olivier BG, Bergmann FT, Ataman M, Babaei P, Bartell JA, Blank LM, Chauhan S, Correia K, et al. 2020. MEMOTE for standardized genome-scale metabolic model testing. Nat Biotechnol 38:272–276. https://doi.org/10.1038/s41587-020-0446-y

67. Ebrahim A, Lerman JA, Palsson BO, Hyduke DR. 2013. COBRApy: Constraints-Based Reconstruction and Analysis for Python. BMC Syst Biol 7:74. https://doi.org/10.1186/1752-0509-7-74

68. Rozas J, Ferrer-Mata A, Sánchez-DelBarrio JC, Guirao-Rico S, Librado P, Ramos-Onsins SE, Sánchez-Gracia A. 2017. DnaSP 6: DNA sequence polymorphism analysis of large data sets. Mol Biol Evol 34:3299–3302. https://doi.org/10.1093/molbev/msx248

69. Dixon P. 2003. VEGAN, a package of R functions for community ecology. J Vegetat Sci 14:927–930. https://doi.org/10.1111/j.1654-1103.2003.tb02228.x

70. Coordinators NR. 2016. Database resources of the National Center for Biotechnology Information. Nucleic Acids Res 44:D7–D19. https://doi.org/10.1093/nar/gkv1290

71. Wood DE, Lu J, Langmead B. 2019. Improved metagenomic analysis with Kraken 2. Genome Biol 20:257. https://doi.org/10.1186/s13059-019-1891-0

72. Price DRG, Bartley K, Blake DP, Karp-Tatham E, Nunn F, Burgess STG, Nisbet AJ. 2021. A Rickettsiella endosymbiont is a potential source of essential B-vitamins for the poultry red mite, *Dermanyssus gallinaceus*. Front Microbiol 12:695346. https://doi.org/10.3389/fmicb.2021.695346

73. Azagi T, Dirks RP, Yebra-Pimentel ES, Schaap PJ, Koehorst JJ, Esser HJ, Sprong H. 2022. Assembly and comparison of Ca. Neoehrlichia mikurensis genomes. Microorganisms 10:1134. https://doi.org/10.3390/microorganisms10061134

74. Garcia-Vozmediano A, Tomassone L, Fonville M, Bertolotti L, Heylen D, Fabri ND, Medlock JM, Nijhof AM, Hansford KM, Sprong H, Krawczyk AI. 2022. The genetic diversity of *Rickettsiella* symbionts in *Ixodes ricinus* throughout Europe. Microb Ecol 84:613–626. https://doi.org/10.1007/s02482-021-01869-7

75. Sassera D, Lo N, Epis S, D'Auria G, Montagna M, Comandatore F, Horner D, Peretó J, Luciano AM, Franciosi F, Ferri E, Crotti E, Bazzocchi C, Daffonchio D, Sacchi L, Moya A, Latorre A, Bandi C. 2011. Phylogenomic evidence for the presence of a flagellum and *cbb3* oxidase in the free-living mitochondrial ancestor. Mol Biol Evol 28:3285–3296. https://doi.org/10.1093/molbev/msr159

76. Cronan JE, Thomas J. 2009. Bacterial fatty acid synthesis and its relationships with polyketide synthetic pathways. Methods Enzymol 459:395–433. https://doi.org/10.1016/S0076-6879(09)04617-5

77. Jain C, Rodriguez-R LM, Phillippy AM, Konstantinidis KT, Aluru S. 2017. High-throughput ANI analysis of 90K prokaryotic genomes reveals clear species boundaries. bioRxiv. https://doi.org/10.1101/225342

78. Diop A, El Karkouri K, Raoult D, Fournier P-E. 2020. Genome sequence-based criteria for demarcation and definition of species in the genus *Rickettsia*. Int J Syst Evol Microbiol 70:1738–1750. https://doi.org/10.1099/ijsem.0.003963

79. Helminiak L, Mishra S, Kim HK. 2022. Pathogenicity and virulence of *Rickettsia*. Virulence 13:1752–1771. https://doi.org/10.1080/21505594.2022.2132047

80. Reed SCO, Lamason RL, Risca VI, Abernathy E, Welch MD. 2014. *Rickettsia* actin-based motility occurs in distinct phases mediated by different actin nucleators. Curr Biol 24:98–103. https://doi.org/10.1016/j.cub.2013.11.025

81. Pillonel T, Bertelli C, Aeby S, de Barsy M, Jacquier N, Kebbi-Beghdadi C, Mueller L, Vouga M, Greub G. 2019. Sequencing the obligate intracellular *Rhabdochlamydia helvetica* within its host *Ixodes ricinus* to investigate their symbiotic relationship. Genome Biol Evol 11:1334–1344. https://doi.org/10.1093/gbe/evz072

82. Guizzo MG, Hatalová T, Frantová H, Zurek L, Kopáček P, Perner J. 2023. *Ixodes ricinus* ticks have a functional association with *Midichloria*
*mitochondrii*. Front Cell Infect Microbiol 12. https://doi.org/10.3389/fci
mb.2022.1081666

83. Cafiso A, Sassera D, Romeo C, Serra V, Hervet C, Bandi C, Plantard O, Bazzocchi C. 2019. *Midichloria mitochondrii*, endosymbiont of *Ixodes ricinus*: evidence for the transmission to the vertebrate host during the tick blood meal. Ticks Tick Borne Dis 10:5–12. https://doi.org/10.1016/j.ttbdis.2018.08.008

84. Gofton AW, Oskam CL, Lo N, Beninati T, Wei H, McCarl V, Murray DC, Paparini A, Greay TL, Holmes AJ, Bunce M, Ryan U, Irwin P. 2015. Inhibition of the endosymbiont *Candidatus Midichloria mitochondrii*" during 16S rRNA gene profiling reveals potential pathogens in *Ixodes* ticks from Australia. Parasit Vectors 8:345. https://doi.org/10.1186/s13071-015-0958-3

85. Ninio C, Plantard O, Serra V, Pollera C, Ferrari N, Cafiso A, Sassera D, Bazzocchi C. 2015. Antibiotic treatment of the hard tick *Ixodes ricinus*: influence on *Midichloria mitochondrii* load following blood meal. Ticks Tick Borne Dis 6:653–657. https://doi.org/10.1016/j.ttbdis.2015.05.011

86. Mariconti M, Epis S, Sacchi L, Biggiogera M, Sassera D, Genchi M, Alberti E, Montagna M, Bandi C, Bazzocchi C. 2012. A study on the presence of flagella in the order Rickettsiales: the case of "Candidatus Midichloria mitochondrii". Microbiology (Reading) 158:1677–1683. https://doi.org/10.1099/mic.0.057174-0

87. Coipan EC, Jahfari S, Fonville M, Maassen CB, van der Giessen J, Takken W, Takumi K, Sprong H. 2013. Spatiotemporal dynamics of emerging pathogens in questing *Ixodes ricinus*. Front Cell Infect Microbiol 3:36. https://doi.org/10.3389/fcimb.2013.00036

88. Kagemann J, Clay K. 2013. Effects of infection by *Arsenophonus* and *Rickettsia* bacteria on the locomotive ability of the ticks *Amblyomma americanum*, *Dermacentor variabilis*, and *Ixodes scapularis*. J Med Entomol 50:155–162. https://doi.org/10.1603/me12086

89. Maitre A, Wu-Chuang A, Mateos-Hernández L, Foucault-Simonin A, Moutailler S, Paoli J-C, Falchi A, Díaz-Sánchez AA, Banović P, Obregón D, Cabezas-Cruz A. 2022. *Rickettsia helvetica* infection is associated with microbiome modulation in *Ixodes ricinus* collected from humans in Serbia. Sci Rep 12:11464. https://doi.org/10.1038/s41598-022-15681-x

90. Haglund CM, Choe JE, Skau CT, Kovar DR, Welch MD. 2010. Rickettsia Sca2 is a bacterial formin-like mediator of actin-based motility. Nat Cell Biol 12:1057–1063. https://doi.org/10.1038/ncb2109

91. Oliver JD, Burkhardt NY, Felsheim RF, Kurtti TJ, Munderloh UG. 2014. Motility characteristics are altered for *Rickettsia bellii* transformed to overexpress a heterologous *rickA* gene. Appl Environ Microbiol 80:1170–1176. https://doi.org/10.1128/AEM.03352-13

92. Sallay B, Vaculová T, Derdáková M, Rusňáková Tarageľová V, Špitalská E, Škultéty L. 2017. Two mice models for transferability of zoonotic bacteria via tick vector. Acta Virol 61:372–376. https://doi.org/10.4149/a/v_2017_319

93. Burgdorfer W, Aeschlimann A, Peter O, Hayes SF, Philip RN. 1979. *Ixodes ricinus*: vector of a hitherto undescribed spotted fever group agent in Switzerland. Acta Trop 36:357–367.

94. Watanabe K, Yukuhiro F, Matsuura Y, Fukatsu T, Noda H. 2014. Intrasperm vertical symbiont transmission. Proc Natl Acad Sci USA 111:7433–7437. https://doi.org/10.1073/pnas.1402476111

95. Hurst LD. 1990. Parasite diversity and the evolution of diploidy, multicellularity and anisogamy. J Theor Biol 144:429–443. https://doi.org/10.1016/s0022-5193(05)80085-2

96. Wang M, Zhu D, Dai J, Zhong Z, Zhang Y, Wang J. 2018. Tissue localization and variation of major symbionts in *Haemaphysalis longicornis*, *Rhipicephalus haemaphysaloides*, and *Dermacentor silvarum* in China. Appl Environ Microbiol 84:e00029-18. https://doi.org/10.1128/AEM.00029-18

97. Hayes SF, Burgdorfer W, Aeschlimann A. 1980. Sexual transmission of spotted fever group rickettsiae by infected male ticks: detection of rickettsiae in immature spermatozoa of *Ixodes ricinus*. Infect Immun 27:638–642. https://doi.org/10.1128/iai.27.2.638-642.1980

98. Gu C, Kim GB, Kim WJ, Kim HU, Lee SY. 2019. Current status and applications of genome-scale metabolic models. Genome Biol 20:121. https://doi.org/10.1186/s13059-019-1730-3

99. Andersson SG, Zomorodipour A, Andersson JO, Sicheritz-Pontén T, Alsmark UC, Podowski RM, Näslund AK, Eriksson AS, Winkler HH, Kurland CG. 1998. The genome sequence of *Rickettsia prowazekii* and the origin of mitochondria. Nature 396:133–140. https://doi.org/10.1038/24094

100. Walker DH, Yu X-J. 2005. Progress in rickettsial genome analysis from pioneering of *Rickettsia prowazekii* to the recent *Rickettsia typhi*. Ann N Y Acad Sci 1063:13–25. https://doi.org/10.1196/annals.1355.003

101. Renesto P, Ogata H, Audic S, Claverie J-M, Raoult D. 2005. Some lessons from *Rickettsia* genomics. FEMS Microbiol Rev 29:99–117. https://doi.org/10.1016/j.femsre.2004.09.002

102. Buysse M, Floriano AM, Gottlieb Y, Nardi T, Comandatore F, Olivieri E, Giannetto A, Palomar AM, Makepeace BL, Bazzocchi C, Cafiso A, Sassera D, Duron O. 2021. A dual endosymbiosis supports nutritional adaptation to hematophagy in the invasive tick *Hyalomma marginatum*. *Elife* 10:e72747. https://doi.org/10.7554/eLife.72747