Molecular characterization of Cardinium, Rickettsia, Spiroplasma and Wolbachia in mite species from citrus orchards

Tetranychidae spider mites are considered key citrus pests in some production areas, especially Tetranychus urticae Koch. Over the past decades, pesticide overuse seems to have promoted T. urticae population selection in citrus orchards. However, the microbiota has also been pointed out as a plausible explanation for population structure or plant host specialisation observed in several arthropod species. In this work, we have determined the incidence of Cardinium, Rickettsia, Spiroplasma and Wolbachia as representatives of major distorter bacteria genera in Aplonobia histricina (Berlese), Eutetranychus banksi (McGregor), Eutetranychus orientalis (Klein), Panonychus citri (McGregor), Tetranychus evansi Baker and Pritchard, Tetranychus turkestani Ugarov and Nikolskii, and T. urticae populations from Spanish citrus orchards. Only Wolbachia was detected by PCR. The multilocus alignment approach and phylogenetic inference indicated that all detected Wolbachia belong to supergroup B. The deep analysis of each 16S rDNA, ftsZ and wsp gene sequences allowed identifying several phylogenetically different Wolbachia sequences. It probably indicates the presence of several different races or strains, all of them belonging to supergroup B. The wsp sequence typing analysis unveiled the presence of the two already identified alleles (61 and 370) and allowed to contribute with five new alleles, supporting the presence of different but related B-races in the studied mite populations. The results are discussed and related to T. urticae population structure, previously observed in Spanish citrus orchards.


Introduction
Spider mites of the family Tetranychidae comprise more than 1300 phytophagous species, out of which about 10% are considered agricultural pests and approximately 10 are key pests of economically important crops (Migeon and Dorkeld 2019). Tetranychus urticae Koch is the most widespread mite, considered one of the citrus key pests, together with the Mediterranean medfly Ceratitis capitata (Wiedemann) and the diaspidid scale Aonidiella aurantii (Maskell) . Phytoseiidae mites-either naturally present in the tree canopy or/and ground cover, or introduced-are the natural enemies providing the biological control of these Tetranychidae, which in integrated pest management can be complemented with a rational application of pesticides (Jacas et al. 2010). However, due to past abuse of pesticides, T. urticae populations in Mediterranean citrus orchards have shown a genetic structuring, which could be attributed to pesticide-driven selection (Pascual-Ruiz et al. 2014). Besides, other studies indicated the presence of selective mating forces or maternal factors that link T. urticae populations' genetic structure to plant host species, which could explain the genetic structuring observed (Marinosci et al. 2015;Aguilar-Fenollosa et al. 2016;Sato et al. 2016). Such forces/factors have remained unsolved in citrus mites of Spain.
By the mid-1960s, bacterial and yeast symbionts of arthropods and nematodes were highlighted as maternal factors affecting the ecology, evolution and reproductive biology of their hosts (Buchner 1965). Over the past 2 decades, this microbiota has become the focus of numerous studies, going from an ecological to a genomic perspective. More recently, the outcomes of these studies are being devised as a new form of biological control, by inducing reproductive barriers with the natural populations mediated by bacterial species (Zabalou et al. 2004;Atyame et al. 2011;Zhou and Li 2016). For example, cytoplasmic incompatibility (CI), a reproductive modification caused by some bacteria, can be used as a population suppression strategy, analogous to the sterile insect technique (SIT) that reduces or eliminates the population, or/and as population replacement, using the bacteria as a vehicle to drive desired phenotypes into natural populations (Brelsfoard and Dobson 2009). Cardinium, Rickettsia, Spiroplasma and Wolbachia are the representative genera of these bacterial distorters that infect many arthropod species (Jeyaprakash and Hoy 2000;Engelstadter and Hurst 2009;Duron and Hurst 2013). Cardinium encompasses a bacterial genus of Bacteroidetes that induces reproductive alterations in its hosts such as CI, parthenogenesis and feminisation Gotoh et al. 2007a;Zhu et al. 2012). Rickettsia and Wolbachia genera belong to Rickettsiales (within alpha-proteobacteria), forming two isolated clades that also induce reproductive alterations (as male feminisation, thelytokous parthenogenesis, CI and male death) and have also been related to pesticide resistance development (Werren 1997;Stouthamer et al. 1999;Stevens et al. 2001;Perlman et al. 2006;Hosokawa et al. 2010;Liu and Guo 2019). Spiroplasma belongs to the Mollicutes (within Firmicutes) and is also involved in the protection of its host against biotic and abiotic stresses (Bolanos et al. 2015;Heyworth and Ferrari 2015;Frago et al. 2017;Guidolin et al. 2018). Recent estimations of arthropod bacterial infestation reached up to 13% for Cardinium, 24% for Rickettsia, 5-10% for Spiroplasma and to 52% for Wolbachia (Duron et al. 2008;Weinert et al. 2015;Mathé-Hubert et al. 2019).
These four genera are transmitted mainly vertically, from mother to offspring, by transovarial infection of eggs. Horizontal transfer has also been reported, either plantmediated or transmitted by some parasitoid species (Russell et al. 2003;Sintupachee 1 3 et al. 2006;Oliver et al. 2010;Ahmed et al. 2015;Li et al. 2017). Due to their intracellular lifestyle (except for some Spiroplasma species), most of these bacteria cannot be grown outside their arthropod host and their identification depends on the application of molecular methods. Whereas bacterial species' identification relies on the positive amplification with species-specific primers, located mainly in the multicopy 16S rDNA locus, the Wolbachia incompatibility strain assignment is performed by multipleloci sequence alignment analysis (MLSA) and phylogenetic inference against reference strains (Russell et al. 2003;Ros et al. 2009). To date, 16 Wolbachia supergroups (named with letters from A to Q, with some recombination events) have been established based on these MLSA analyses (Lo et al. 2002(Lo et al. , 2007Bordenstein and Rosengaus 2005;Ros et al. 2009;Augustinos et al. 2011;Pascar and Chandler 2018).
As indicated previously, some of these bacterial species are involved in CI (being able to modulate population genetic structure), pesticide resistance and biotic/abiotic stress resistance (water and temperature). Therefore, the determination of their presence in the natural populations of Tetranychidae is important to ascertain how they may affect the host population structure.
In this work, we studied the incidence and frequency of infection of  Table 1 lists the mites (mainly Tetranychidae) collected mainly from Spanish citrus orchards or from laboratory rearing colonies, and the insect species used as positive controls for PCR. The numbers of specimens per species or population are also included in Table 1.

DNA extraction and verification
Total DNA was extracted from isolated, ethanol-washed specimens following a modified 'salting out' protocol (Pérez-Sayas et al. 2015). Briefly, each surface-disinfected specimen was air-dried, isolated in a 1.5-ml Eppendorf tube and crushed in TNS + Prot-K solution at 60 °C; proteins were precipitated with 5 M NaCl by centrifugation and the nucleic acid fraction was precipitated with 2-propanol. The extracted DNA from non-Acari specimens was quantified with Nanodrop 2000 (Thermo Scientific, Wilmington, DE, USA). The Acari specimens' DNA extractions were subjected to PCR with 18SrDNA primers (see Table 2) to ascertain the presence of DNA, as previously done with minute specimens (Pérez-Sayas et al. 2015). ) and references, used to determine the incidence of bacterial symbionts in our samples Target Primer name

Cardinium, Rickettsia, Spiroplasma and Wolbachia diagnostic PCR
The incidence of each bacterial symbiont was determined by positive PCR reactions with specific primers (listed in Table 2), targeting the 16S rDNA in each specimen collected. Due to the limiting factor of Acari source DNA, a secondary specific amplification was devised over a first (primary) amplification of whole 16S rDNA fragment, using the universal primers listed in Table 2, as devised for other insect-bacteria groups (van Ham et al. 1997;Russell et al. 2003). The primary PCR was performed using 1 µl of DNA extraction, whereas the specific secondary and diagnostic PCR was performed with 1-2 µl of the primary PCR. Amplification conditions varied slightly between bacterial species (see Table 2 for reaction volume, magnesium concentration and annealing temperatures), using 1 U of FIREPol polymerase (Solis BioDyne, Tartu, Estonia) with the appropriate 1× buffer, with 0.2 mM dNTPs and 0.4 mM of each primer. Amplification was performed in a C1000 Bio-Rad thermocycler (Applied Biosystems, Foster City, CA, USA) under the following amplification conditions: a first denaturing step at 92-95 ºC for 2-5 min, followed by 30-40 cycles of 92-95 ºC for 30 s, 52-58 ºC for 30 s and 72 ºC for 30-60 s, with a final extension at 72 ºC for 5 min (see Supplementary information). For each amplification run, at least one negative control (ultrapure water added instead of DNA sample) and one positive specimen (of the species listed in Table 1; at least one per symbiont species to be determined) were included to ascertain the false positives (either due to contaminated reagents or environmental contamination) and negatives (due to failure of amplification or low DNA concentration), respectively. Amplification was verified by agarose gel (2% low EEO DA Agarose, Pronadisa, Sumilab, Madrid, Spain) electrophoresis in 1× TAE, stained with Gel-Red (Biotium, Hayward, CA, USA). Single, expected-size PCR fragments were considered positive when matching the size of the positive controls. Each specimen was considered harbouring Cardinium, Rickettsia, Spiroplasma or Wolbachia, when at least two PCR reactions give positive results of the three performed.
Positive PCR fragments were independently purified with Illustra ExoStar (GE Healthcare Life Sciences, Chalfont St. Giles, UK) following the manufacturer's recommendations. Bidirectional Sanger sequencing using Bigdye terminator v.3.1 cycle sequencing kit (Thermo Fisher Scientific, Vilnius, Lithuania) with each amplification primer was performed at the Sequencing service of the University of Valencia (Servei Central de Suport a la Investigació Experimental [SCSIE], Universitat de València, Spain), following the manufacturer's instructions. Reactions were run in an ABI 3730XL DNA analyser (Thermo Fisher Scientific, Carlsbad, CA, USA) following the manufacturer's instructions.

Wolbachia wsp and ftsZ amplification and sequencing
To assign Wolbachia into the established supergroups, we used the MLSA approach by amplifying and sequencing the genes corresponding to cell division protein FtsZ (ftsZ) and the Wolbachia surface protein (wsp), in addition to the 16S rDNA described above (Braig et al. 1998;Zhou et al. 1998;Lo et al. 2002;Casiraghi et al. 2005;Baldo et al. 2006). PCRs were conducted independently using 1-2 µl of undiluted specimen DNAs with primers and conditions, as listed in Table 2, using 1 U of FIREPol polymerase with 0.2 mM dNTPs and 0.4 mM of each primer for the 16S rDNA amplification. Amplifications were performed in a Bio-Rad thermal cycler with the following amplification conditions: a first denaturing step at 94-95 ºC for 2-5 min, followed by 36-40 cycles of 94-95 ºC for 30 s, 54-55 ºC for 1 3 45-60 s, and 72 ºC for 60-90 s, with a final extension at 72 ºC for 5 min (see Supplementary information). Similarly, positive (other arthropods specimens harbouring known types of Wolbachia and/or Wolbachia positive T. urticae samples) and negative (DNA-free PCR mixture) controls were included in each amplification run. Positive PCR fragments were purified as described above and sequenced bidirectionally with amplification primers, at the same SCSIE sequencing service.

Sequence analysis
The consensus sequence for each PCR product was obtained using the programme STADEN Package (Staden 1996). Consensus sequences were blasted against the nonredundant database to confirm fragment identity prior to alignment construction (BLAST; Altschul et al. 1997).
16S ribosomal DNA, ftsZ and wsp obtained consensus sequences and those retrieved from databases were independently aligned using CLUSTALW (as in MEGA X; Kumar et al. 2018) (for 16S rDNA) or with GENEDOC (Nicholas andNicholas 1994-1998). In GENEDOC, we used Blosum62 score table for coding regions ftsZ and wsp, whereas for 16S rDNA we used PAM 65 score table, setting alignment cost at 20 for constant length, 8 for gap opening and 4 for gap extension (for ftsZ and wsp, alignment was performed with translated sequences, re-gapping the nucleotide alignments). Moreover, 16S rDNA, ftsZ and wsp consensus sequences were concatenated in a single FASTA file previously to perform the multilocus sequence alignment (MLSA). Outgroups were retrieved from the databases and sequences corresponding to the same species were concatenated in the same order as the MLSA (Table S1).
The wsp sequences were assigned to the corresponding allelic profile by comparing the four hypervariable regions (HVRs) against the Wolbachia wsp multilocus sequence typing (MLST) database (https ://pubml st.org/wolba chia/ [last accessed 10/March/2020]; Baldo et al. 2005). Novel allele sequences were submitted to the database curators for their inclusion as new alleles after they registered as new sequences in NCBI.
Gene tree inference was conducted in MEGA X, after determining the best-fit evolutionary distance model (GTR) for each gene alignment and for the MLSA, as implemented in MEGA X. Bayesian phylogenies were obtained using a Markov Chain Monte Carlo (MCMC) method implemented in BEAST v.1.10.4 programme . BEAST output was analysed using TRACER v.1.7.1, applying values of more than 200 of the effective sample size (ESS) (Rambaut et al. 2018). A maximum clade credibility tree was generated after burning 10% samples with posterior probability limit > 0.5 using TreeAnnotator, as implemented in BEAST. Species phylogroups were defined by a posterior probability > 0.95 using referenced strains, known to belong to these groups.

Data availability
All new sequences have been deposited in GenBank from MN123012 to MN123230 for 16S rDNA, MN187577-MN187703 for wsp gene region, and MN187704-MN187866 for ftsZ gene region (see Table S1 for the complete list).

Incidence of Cardinium, Rickettsia, Spiroplasma and Wolbachia in mite populations
Cardinium, Rickettsia and Spiroplasma species-specific primer pairs gave negative results in all the mite species and populations tested, despite their amplification efficacy being positive with the corresponding arthropod control samples. All samples were tested for DNA presence, as routinely done with such minute specimens, by amplification of 18S rDNA (Pérez-Sayas et al. 2015). Only the 16S universal and Wolbachia-specific primers (either 16S rDNA, ftsZ or wsp) rendered positive results. Wolbachia was present in almost all the mite species and/or populations tested with a prevalence ranging from 10 to 100% (Fig. 1), as previously reported (Zug and Hammerstein 2012;Weinert et al. 2015;Zhu et al. 2018). The exception was P. citri, which showed a prevalence of 0-10%. This is the first time that Wolbachia is reported in this mite species (Zélé et al. 2018a, b;Zhu et al. 2018).
Other authors have detected double infections in Tetranychus species; two of them were included in our study, namely T. urticae and T. evansi (Enigl and Schausberger 2007;Weinert et al. 2009;Xie et al. 2016;Staudacher et al. 2017;Zélé et al. 2018a, b). These studies found that Tetranychus truncatus Ehara showed the combinations Wolbachia and Cardinium or Spiroplasma and Rickettsia, whereas T. evansi, Tetranychus ludeni Zacher and T. urticae showed only the Wolbachia and Cardinium combinations (Zhang et al. , 2016Zélé et al. 2018a). Indeed, the double infection Wolbachia and Cardinium (W + C) is the most common in Tetranychidae (Zélé et al. 2018a, b). All these studies used two specific primer pairs: one pair targeting the 16S rDNA in Cardinium and Spiroplasma and the second one targeting a species-specific gene (i.e., gyrB for Cardinium, rpoB for Note that only one population of Panonychus citri showed Wolbachia infection (0.87%). The total number of individuals tested are, in order of species appearance and from left to right: 10, 9, 15, 115, 80, 36, 26 and 368. Spiroplasma and gltA for Rickettsia). As indicated, the diagnostic primers used, except the Cardinium ones, differ from other mite working groups but render positive results with arthropod species used as positive controls (see Table 2 for references of each primer pair).
Our aim was to detect each bacterial species based on the same target gene to include all of them in a phylogenetic study to determine the presence of more than one strain, for which our primer selection based on previous works targeting the multicopy gene 16S rDNA. In previous studies, we have observed a population structure in T. urticae within Spanish populations (by microsatellite analyses), that may be attributed to different Wolbachia operational taxonomic units, as noticed in the present study (operational taxonomic units as described for 16S rDNA sequences diverging more than 3-5% as in microbiome analyses; see below) (Aguilar-Fenollosa et al. 2012;Pascual-Ruiz et al. 2014). However, as Wolbachia was the only bacterial reproductive distorter detected in our study, a deep analysis of the Wolbachia sequences obtained was required to clarify the situation.

Phylogeny of Wolbachia and strain identification
Wolbachia is a group of bacterial strains that can be assigned to supergroups following a multilocus phylogeny approach, as indicated previously (Ros et al. 2009). The sequence of coding genes of the cell division protein FtsZ (ftsZ gene) or the Wolbachia surface protein (wsp) are routinely used for placement of Wolbachia strains into the established supergroups A to K Gotoh et al. 2003;Casiraghi et al. 2005;Baldo et al. 2007), whereas the 16S rDNA, is routinely used to determine bacterial species identity in microbiome studies.
Here, we have estimated the tree phylogeny of Wolbachia from several mite species with a multilocus alignment (MLSA) of concatenated 16S rDNA, wsp and ftsZ, and then analysed each locus independently by different tree-reconstruction methods (Maximum Parsimony [MP], Maximum Likelihood [ML] and Bayesian). Using the MLSA approach, either by MP (MP was used to compare against precedent work by Zhang et al. 2013), ML or by Bayesian inference, almost all the mites' new Wolbachia sequences clustered within the supergroup B, except the ones from the Brazilian population of T. evansi (TeBr45 and TeBr70) that clustered either basal in the B group (Bayesian inference; Fig. 2a) or between A-, K-, C-Wolbachia supergroups (ML) (Fig. 2b).
Due to the scarcity of DNA material obtained from these minute mites, it was impossible to obtain a sequencing grade wsp fragment for some samples, which limited the number of samples used for this MLSA to 90 individuals. Consequently, the power of MLSA to determine the presence of more than one Wolbachia strain in our samples was limited. When we analysed the MLSA solely composed by 16S rDNA + ftsZ fragments, we increased the figure to 121 newly concatenated sequences, despite that for the tree inference 100% identical sequences from the same mite population were removed to reduce the computing time (Fig. 3). Limiting data and samples reduced the resolution of the trees and improved deep branching in some cases, whereas in others, and due to positive selection Fig. 2 Phylogenetic inference from MLSA (containing 16S rDNA, ftsZ and wsp, concatenated sequences) of 72 Wolbachia specimens (indicated with the corresponding species name, sample code and GenBank accession number of the 16S sequence of each specimen) inferred using the following method: a Bayesian analysis or b Maximum Likelihood (ML). Evolutionary distances were computed using the Tamura 3-parameter method after modelling rate variation among sites with a gamma distribution (shape parameter = 0.08). All ambiguous positions were removed for each sequence pair. Evolutionary analyses were conducted in MEGA X. Wolbachia supergroups are indicated in the krone section outside with different patterns.
▸ detected in some wsp lineages, clustering of sequences belonging to the same supergroups did not match previous works (Schulenburg et al. 2000;Ros et al. 2009). In this case (16S rDNA + ftsZ), all Tetranychidae sequences clustered together within the supergroup B, not showing any structuring between the Brazilian population of T. evansi nor the already characterised as different members of supergroup B (including in this last group all the B-Wolbachia from various insect species with different reproductive modes) (Ros et al. 2009(Ros et al. , 2012. Further, when each gene fragment was independently analysed, we could observe a supported differentiation that depends on the fragment type (coding or non-coding). Our limiting sequence (by the number of samples and available supergroups), wsp, gave different tree inferences ( Fig. 4 and Fig. S1), keeping in both cases supergroups A, C and E as basal with high posterior probabilities or bootstrap values. While the B-supergroup was split into three clusters (B1, B2 and B3; Fig. 4), the first two, B1 and B2, included many of the outgroup sequences. Some of them were linked but not completely isolating species with feminisation or thelytoky reproductive specialisations. Group B3 included many species (from outgroup) with identified CI, with all of our sequences (van Meer et al. 1999). Despite this, group B3 seemed to also show an internal split into three other groups with posterior probabilities higher than 0.96; the results did not find any relationship between Wolbachia taxonomic unit (B-sub-sub-strain; sequences that show high-sequence divergence, conforming differential taxonomic units) and host plant or mite populations, as previously found with microsatellites. Further, ftsZ phylogenies placed A-supergroup sequences in a basal cluster to B-supergroup, which is subdivided into three subgroups (B1-B3 in Fig. 5 and Fig. S2), on which again sequences of T. evansi from Brazil roots in the most basal subgroup (B1). The 16S rDNA phylogenies were most resolute, supporting the clustering of supergroups, as previously published ( Fig. 6 and Fig. S3) (Gotoh et al. 2003(Gotoh et al. , 2007bRos et al. 2009Ros et al. , 2012Suh et al. 2015). With this marker, supergroup B was split into five subgroups (B1-B5 in Fig. 6), with T. evansi Brazil population sequences mostly concentrated within subgroup B3. In this phylogenetic reconstruction, B-Wolbachia from vector insects like Bactericera cockerelli (Šulc) (Hemiptera: Triozidae) (EF372596) and Diaphorina citri Kuwayama (Hemiptera: Psyllidae) (GU563892) or other pests like Naupactus cervinus Boheman (Coleoptera: Curculionidae) (GQ402143) or mites like Bryobia spp. (i.e., EU499318) were clustered together in a well-supported clade B2. However, the T. urticae T2 reference sequence (EU499319) clustered within subgroup B5, which contained some populations of T. urticae, including those from our previous studies on which a genetic structure was devised (Aguilar-Fenollosa et al. 2012). Group B3 contained samples of T. turkestani and the majority of T. evansi Brazil population sequences. The sequence divergence of 16S rDNA among these subgroups was sometimes higher than the reference 3% used in microbiome analysis, indicating that this clustering reflects the diversity of Wolbachia races within Tetranychidae mites .
In addition to these phylogeny-based classification methods (MLSA or 16S rDNA barcoding), other methods to identify Wolbachia strains have been developed in other studies Fig. 3 Phylogenetic inference from MLSA (containing only 16S rDNA and ftsZ concatenated sequences) of 90 Wolbachia specimens (indicated with the corresponding species name, sample code and GenBank accession number of the 16S sequence of each specimen) was inferred using a Bayesian analysis or b Maximum Likelihood (ML). Evolutionary distances were computed using the Tamura 3-parameter method after modelling rate variation among sites with a gamma distribution (shape parameter = 0.27). Evolutionary analyses were conducted in MEGA X. Wolbachia supergroups are indicated in the krone section outside with different patterns.
1 3 (revised in Bleidorn and Gerth 2018). One of them is the MLST system, based on allele assignment of gatB, coxA, hcpA, fbpA and ftsZ genes (allele assignment was per single nucleotide difference with reference strain in a concatenated sequence of these five genes) (Baldo et al. 2006;Jolley and Maiden 2010). As we only sequenced gene ftsZ, we could not use the whole MLST approach; however, based on this kind of study, all T. urticae specimens (ours and some already characterised as different) were assigned to the ftsZ locus 23. Whereas the Brazilian population of T. evansi presented the ftsZ locus 179. Recently, the same authors included the allele typing with only wsp gene due to its key features (single-copy gene, present in all Rickettsiales order, with evidence of strong stabilising selection and generally used as phylogenetic marker) and matching one of our sequenced genes (Baldo et al. 2006;Jolley and Maiden 2010). Following this wsp sequence typing, we were able to assign our B-Wolbachia sequences to different wsp alleles, including the description of five new wsp loci (submitted to the MLST database on 9 March 2020, three presented here as X1 to X3). The wsp locus 61 (HVR1:18; HVR2:16; HVR3:23; HVR4:16) was the predominant one in almost all T. urticae feeding in citrus (54%; 66 out of 122), followed by wsp locus X1 (24%, HVR1:18; HVR2:16; HVR3:23; HVR4:274) in samples from Festuca arundinacea cover and other populations (24%). Tetranychus urticae feeding . Due to the reduced number of individuals per population tested, we were not able to conduct a proper analysis of diversity. However, we were able to clearly identify different alleles, indicating that there exists more than one strain of Wolbachia in some of our populations.
Considering phylogenies and wsp MLST, we can conclude that T. urticae populations show different B-Wolbachia strains. Their involvement in mite reproduction could explain the T. urticae population structure previously established in Spanish citrus orchards, deserving further research to determine the link between each strain and reproductive isolation (Aguilar-Fenollosa et al. 2012Zhang et al. 2013;Pascual-Ruiz et al. 2014). This result is in line with other studies in which D. citri, one of the vectors of Huanglongbing (HLB), seems to be infected by two B-Wolbachia races, affecting their population structure and differential transmission of Candidatus Liberobacter, the plant pathogenic bacterium causing HLB (Chu et al. 2019). Similarly, T. urticae populations from Korean greenhouses have been reported to harbour two Wolbachia races based on their wsp sequences, showing diverse patterns of CI that matched the host plant as the main phenotypic effect, similar to the population structuration previously devised due to CI in Chinese and Japanese T. urticae populations or in recent invasive events in Europe (Gotoh et al. 2007b;Boubou et al. 2011Boubou et al. , 2012Xie et al. 2011;Zhang et al. 2013;Suh et al. 2015). However, with the samples analysed, we could not relate each identified B-Wolbachia strain (or race) with a specific genome structuration, derived either by pest management, host plant specificity or even by its reproductive alteration pattern, which deserves further study.

Final remarks
We have identified only one bacterial species, Wolbachia, of the four manipulative tenant bacteria tested in our mite target populations. This bacterial species was assigned by phylogenetic analysis to the B-supergroup, highlighting the existence of several races or strains within them. Sequence typing of wsp gene allowed the assignment to several alleles (mainly alleles 61 and 370) and the description of five new alleles. The presence of several strains could be explained by the biology of Wolbachia, either by an effect in Fig. 6 Phylogenetic inference, using only the 16S region of 126 Wolbachia specimens (indicated with the corresponding species name, sample code and GenBank accession number), was performed using the Bayesian analysis under the GTR + I + Γ model of DNA substitution. Wolbachia supergroups are indicated in the krone section outside with different patterns. the host reproductive strategy (population isolation) or by recent invasive events. Both hypotheses require further study.