Abstract
The continuously high incidence of some endosymbionts in arthropods despite potential conflicts with their hosts is often explained by obligatory relationships, in which the host is fully dependent on its endosymbiont, fitness advantages conferred on hosts by facultative endosymbionts, or reproductive manipulation of hosts by endosymbionts (typically facultative). Yet continuously endosymbiont high incidence is sometimes observed without clear evidence supporting any of these mechanisms. This situation could potentially be explained by the presence of several coinfecting strains of the same endosymbiont species, each affecting the host differently such that their effects counteract one another when studied collectively. Here, we investigated Wolbachia endosymbionts of fleas, which stably persist in high loads in all females, with no indication that any of the above mechanisms explain their continuously high incidence. We sequenced fleas and identified two Wolbachia strains, designated as wSc1 and wSc2. We then correlated the strain composition in fleas with measures of their reproductive success. We found that fleas with high wSc1 and low wSc2 loads had a higher reproductive success than fleas that had high loads of both strains, low loads of both strains, or no Wolbachia, suggesting that wSc1 may provide a direct fitness advantage to their hosts. Conversely, the number of males and total offspring was negatively correlated with wSc2 levels, supporting male killing. Our research demonstrates that endosymbionts’ continuously high incidence may persist through intricate relationships in nature.
Similar content being viewed by others
Introduction
Endosymbiont incidence in arthropod host populations depends on factors such as endosymbiont transmission modes, interactions with the hosts and other endosymbionts, and effects on the host behavior and reproduction, all of which are environmentally influenced1. However, in many cases, despite potential conflicts with their hosts and fluctuating ecological conditions, there exists a continuously high incidence of endosymbionts in arthropod host populations2,3,4. Consistently high incidence is mostly explained by three types of mechanisms5,6. The first involves obligatory relationships, in which the host’s survival, reproduction, or both are entirely dependent on its endosymbiont1,7. The second mechanism can occur even when the host is not fully reliant on the endosymbiont and involves direct fitness benefits. These advantages may include nutrient supplementation or protection from predators or harsh environmental conditions, provided by facultative endosymbionts under certain conditions8,9,10,11,12. The third type of mechanism is reproductive manipulation by the endosymbiont (typically facultative), which enhances its transmission13. Reproductive manipulation may distort the sex ratio of arthropod hosts, for example, when they provide an advantage to female offspring over their non-transmitting male siblings by inducing male killing, parthenogenesis, or feminization. The former phenomenon leads to males’ death, reducing the total offspring number14, and the latter turns genetic males into fertile phenotypic females with functional ovaries15. Another type of reproductive manipulation that does not distort the host sex ratio is cytoplasmic incompatibility (CI), where the offspring of the cross between infected males and uninfected females fail to develop, thereby providing a reproductive advantage to infected females16.
Nevertheless, there are cases of continuously high endosymbiont incidence with unclear drivers (e.g17,18,19,20). One such unusual case is Wolbachia bacteria, which stably persist in high loads in all female Synosternus cleopatrae fleas but not in males, where they occur only sporadically in low loads21. Moreover, despite the endosymbiont being found in all female fleas, we could not find any experimental indications of obligatory relationship, facultative fitness advantage, or reproductive manipulation associated with the endosymbiont6. An alternative explanation for the continuously high incidence of Wolbachia and other endosymbionts in their arthropod hosts with unclear drivers could be that there is coinfection by multiple endosymbiont strains of the same species, a situation that exists in other flea systems22. These strains may have distinct and sometimes opposing effects on their host23, which may become obscured when they are studied collectively. Moreover, it is possible that one of the strains damages male hosts, leading to the sex-bias endosymbiont incidence that is observed in this system.
Two lines of evidence suggest that coinfection can solve at least some of the unexplained cases of continuously high endosymbiont incidence. First, the coinfection of arthropods by multiple endosymbiont strains is common in nature (e.g23,24,25,26,27). Second, recent genomic analyses22,28,29, supported by a handful of experiments30,31,32,33,34,35, provide some evidence for the distinct effects of coinfecting endosymbiont strains on their shared host. These include intriguing experimental examples of Wolbachia strain-specific effects in mosquitoes, moths, flower bugs, wasps, Drosophila, and butterflies30,31,32,33,34,36,37,38,39. Thus, we aimed to determine whether differential effects of coinfecting Wolbachia strains could explain their continuously high Wolbachia incidence in female fleas, despite the lack of experimental support for any of the three proposed mechanisms.
Here we combined field sampling, various molecular techniques, and new analyses of data from a previous experiment that manipulated Wolbachia infection to address the following three, non-mutually exclusive, hypotheses: (H1) S. cleopatrae fleas are coinfected by multiple Wolbachia strains, (H2) these include strains that confer direct fitness advantages, or (H3) cause sex-distortion reproduction manipulations in fleas. We found evidence of coinfections of two Wolbachia strains in both male and female S. cleopatrae, with initial indications that one of these strains confers a direct fitness advantage to the host, while the other induces male killing. These results suggest that contrasting effects of two coinfecting endosymbiont strains may explain how they collectively have a continuously high incidence in host populations.
Materials and methods
Study approach and rationale
The study included multiple steps described in Supplementary Materials, Fig. S1. To test H1, which postulates that the S. cleopatrae fleas are coinfected by multiple Wolbachia strains, we explored Wolbachia genetic variation in S. cleopatrae fleas collected from the field and a laboratory colony (Supplementary Text S1;21). To do this, we subjected the flea DNA extracts to next-generation multilocus sequence typing (NGMLST), targeting Wolbachia’s five housekeeping genes. In most of the fleas, we found two variants per gene. We further assessed whether in these coinfected fleas, there are only two or more—due to possible recombination—Wolbachia variants. We used three strategies to address this question. First, we conducted whole-genome sequencing (WGS) of two flea pools, each including a unique combination of variants across genes but only one variant per gene. Second, on the same flea extracts, we conducted three quantitative real-time PCR (qPCR) assays that targeted the fbpA gene, using either general primers that are designed to amplify all Wolbachia bacteria (hereafter Wolbachia-general) or primers that are specific for each variant. We then compared the Wolbachia-general load and the summed load of the two variant-specific assays. Third, we compared the structures of five phylogenetic trees, each constructed based on the sequences of a single MLST gene. In the final step, we created a phylogenetic tree of concatenated MLST sequences, which allowed us to explore evolutionary relationships between the two Wolbachia variants in S. cleopatrae and in relation to six documented Wolbachia supergroups. Altogether, the results of the above steps suggest that two Wolbachia strains, designated as wSc1 and wSc2, coinfect S. cleopatrae fleas (see Results).
To test H2 and H3 regarding the distinct effects of wSc1 and wSc2 on the fleas, we subjected the flea DNA extracts, which were collected during our previous Wolbachia infection manipulation (Supplementary Text S26;), to strain-specific qPCRs. In short, in the previous study, measures of reproductive success were quantified for four groups of S. cleopatrae. These fleas all originated from the same colony stock but had different combinations of Wolbachia status (Wolbachia-free obtained by antibiotic treatment or Wolbachia-positive) and physiological age (females fed on rodents for either five or ten days6;). More details about the study design of this experiment are provided in the Supplementary Text S2. To test H2 regarding the distinct fitness advantages of the two strains, we correlated their quantities with the fleas’ integrated reproductive success index. Similarly, to test H3 regarding distinct reproductive manipulations of the two strains, we correlated the strain loads with the number of female and male flea offspring (Statistical Analyses). Once we found evidence for reproductive manipulation, we tested whether the mechanism was more likely to be male killing or feminization by correlating the quantities of the two strains and the total number of offspring. A decrease in male offspring with no change in the total number of offspring would indicate feminization40. In contrast, a simultaneous reduction in the number of male offspring and the total number of offspring would indicate male killing41. Note that we used both Wolbachia-free and Wolbachia-positive groups for some of the analyses testing H2, but only the latter group for testing H3.
Sources of field-collected and laboratory fleas
S. cleopatrae fleas were collected during a field survey in 2011 at three sites across the northwestern Negev Desert’s sands in Israel42. In parallel, S. cleopatrae fleas were bred under laboratory conditions on G. andersoni and G. pyramidum rodents6. More details are provided in Supplementary Text S1.
DNA extraction
Fleas used in this study were preserved in 70% ethanol before DNA extraction. DNA was extracted from individual fleas using the QIAGEN DNeasy Blood and Tissue Kit. A negative control extraction was included in each extraction session and all the following analyses.
Next-generation multilocus sequence typing (NGMLST)
We conducted NGMLST analysis on DNA extracts of 128 female fleas: 96 were collected from field-trapped rodents and 32 from the laboratory colony (Supplementary Text S1). We targeted five housekeeping genes commonly used in Wolbachia phylogenetic analysis (Table 1;43). These genes were amplified from each flea sample using PCR and were evenly pooled based on Qubit fluorometer measurements. The resulting 128 pooled samples were subjected to next-generation sequencing. Paired-end 2 × 300 bp reads were generated on the Illumina MiSeq platform at the DNA services facility of the Research Resources Center at the University of Illinois Chicago.
The following bioinformatic analyses were performed on the sequencing results: (i) all primers were removed; (ii) each gene was analyzed—from filtering and trimming to producing a final product amplicon sequence variant table—using the divisive amplicon denoising algorithm (DADA2); (iii) SNPs were considered true if the relative abundance of the relevant sequence was higher than 0.1% within the sample and it occurred in more than 1% of the reads of the specific gene. The bioinformatic analyses were performed in R44 using the dada2 package45.
The coxA and hcpA primers detected less variance than other genes (see Results). To confirm that this was due to a technical limitation and not a biological phenomenon, we redesigned PCR primers based on additional Wolbachia sequences that were derived from the GenBank database (Table 1). The products were amplified with PCR and sequenced with the Sanger method (Macrogen, Europe). The PCR conditions are provided in the Supplementary Text S3.
Quantitative real-time PCRs (qPCRs)
The Wolbachia-general qPCR quantified the total abundance of Wolbachia using the primers and conditions described in Flatau et al. (Table 1;21).
We designed variant-specific (designated as wSc1 and wSc2) qPCRs targeting each MLST gene (Table 1), based on sequences obtained with NGMLST (for gatB, fbpA, and ftsZ) and with PCR (for coxA and hcpA). Standard curves were established from samples with known numbers of wSc1 and wSc2 pUC-GW-Kan plasmids that were synthesized with the concatenated sequences of all five genes from each strain at GENEWIZ from Azenta Life Sciences. They were used to determine the absolute amount of Wolbachia strains within a pooled flea DNA extraction sample.
We confirmed the specificity of the qPCR assays by running all assays on plasmids with only one Wolbachia variant, the other Wolbachia variant, and the pooled flea DNA extraction samples containing both. Even in the highest concentrations, the cross-reaction rate never passed 0.1% for any of the genes except hcpA wSc2, which showed a cross-reaction of 0.6% (Supplementary Table S1).
All specific qPCR assay conditions are provided in the Supplementary Text S3. We found out that all the variant-specific primer pairs distinguish similarly between the two strains. Thus, for the analysis of the samples from the Wolbachia infection manipulation, we used only qPCRs for the detection of wSc1 and wSc2 variants targeting fbpA. Subsequently, to confirm that these two assays covered all Wolbachia bacteria present in tested samples, we compared the cumulative results from both qPCRs on 74 flea samples with the results obtained from the Wolbachia-general qPCR.
Whole genome sequencing of Wolbachia wSc1 and Wolbachia wSc2
We created two pools of flea extracts by assembling all the rare samples with a unique combination of variants across genes, but only one variant per gene. To do this, we gathered hundreds of flea DNA samples in our laboratory collection and subjected them to variant-specific qPCRs for each of the five MLST genes. These resulted in a total of 10 qPCRs per flea extract. After mixing the extracts that had identical MLST profiles into two unique DNA pools, we confirmed their genetic profile by subjecting each pool to another run of variant-specific qPCRs.
Between 50 and 90 ng of gDNA purified from each pool was input into the 2 S Turbo DNA Library Kit (Swift Biosciences). All reactions were carried out at 20% of the manufacturer’s recommended volumes, with dual 6-bp indexes incorporated during the final 12-cycle PCR step. The resulting libraries were pooled and sequenced on an Illumina HiSeq X Ten instrument by Psomagen (Rockville, MD) to generate 151-base paired-end reads. We removed adaptors from the Illumina reads using trimmomatic (v0.39) in paired-end mode with the following settings: four allowed mismatches to the seed, a palindrome clip threshold of 30, a simple clip threshold of ten, and discarding trimmed reads ≤ 30 bases. Paired Illumina reads that survived trimming were mapped to the genomes of Wolbachia bacteria from various arthropods (Wolbachia wCfeJ from Ctenocephalides felis, Wolbachia wDci from Diaphorina citri, and Wolbachia wCon from Tribolium confusum) using bowtie2 (v 2.5.4). These reference genomes were chosen from the pool of Wolbachia whole genomes due to the similarity of their MLST genes to the two strains we found in the S. cleopatrae fleas (Fig. 1). The list of all whole genomes used for comparisons and their NCBI accession numbers are provided in the Supplementary Table S2.
A concatenated phylogenetic tree, based on the five multilocus sequence typing (MLST) genes (gatB, coxA, hcpA, fbpA, and ftsZ) combined into a 1988 bp contig. The tree includes contigs of Wolbachia strains from different hosts (indicated by drawings) and supergroups (capital letters; more details in the Supplementary Table S2). The Wolbachia strains of Synosternus cleopatrae fleas described in this study, are marked in bold.
Phylogenetic trees
Due to the low coverage of the whole genome sequencing data (see Results), we generated maximum-likelihood phylogenetic trees from the NGSMLST sequences. We generated a tree for each MLST target, as well as for concatenated MLST sequences. To create the concatenated tree, we prepared two contigs by concatenating sequences of the five MLST genes from each of the two sets of unique single-strain samples. Similarly, we prepared contigs for all reference Wolbachia bacteria included in the trees.
To cluster the two variants into the established Wolbachia supergroups, we added to the trees reference sequences from six out of the 19 existing supergroups46. These sequences, which represent the most relevant supergroups for arthropod symbionts, namely, supergroups A, B, D, F, and H, were downloaded from the Wolbachia MLST database (PubMLST). We aimed for a maximum of five references per supergroup. When more sequences were available, we maximized the diversity by choosing Wolbachia strains found in different arthropod groups. In addition, to increase the number of reference sequences from supergroups C and D, in which the Wolbachia documented from the C. felis cat flea was originally clustered22, we added the sequences of strains that were documented in nematodes, as well as extracted the Wolbachia MLST sequences from the whole genomes of Wolbachia strains wBp, wLsig, wOo, wOv, and wDimm, which were not available in the PubMLST database. Finally, to include all the other sequences of Wolbachia strains that are associated with fleas, we added the MLST sequences of the other Wolbachia strains documented in C. felis (wCfeT, wCfeJ, and wCfeF) and those reported from C. orientis (wCori) fleas. More details on the isolates used in the tree can be found in Supplementary Table S2.
All trees were created with the Kimura parameters model using the PhyML plugin in Geneious Prime software47. Transition and transversion ratio as well as gamma distribution parameters were estimated, and the proportion of invariable sites was fixed. The number of substitution rate categories was four, and topology, length, and rate were optimized. We aligned the MLST sequences by the Geneious alignment algorithm and constructed maximum likelihood phylogenetic trees with 1000 bootstrap replications.
Statistical analyses
All the statistical analyses were performed using general linear models. To address whether only two Wolbachia variants coinfect the fleas (H1), we correlated the expected overall Wolbachia load (based on Wolbachia-general qPCR) and the sum of loads of the two Wolbachia strains (based on the variant-specific qPCRs).
To address H2 and H3, we performed analyses in which the independent variables were one of two strain-related variables, the maternal physiological age, and the interaction between the strain-related variables and age. Selecting strain-related variables was intended to tackle two distinct inquiries. The associations of the dependent variables with the strain-specific loads were used to test whether the effects are influenced by the load of each strain. The associations of the dependent variables with Wolbachia strain composition would reveal whether the effects are contingent upon a particular composition of the two. For the latter analysis, the strains were assigned a ‘high’ or ‘low’ load relative to the median load of each strain48,49, and the Wolbachia-free group was added.
The dependent variable of the analysis addressing H2 was the integrated index of reproductive success (RS). The index was calculated as follows: \(\:RS={\sum\:}_{i}^{NF}({BS}_{F}\times\:{PS}_{F})+{\sum\:}_{l}^{NM}({BS}_{M}\times\:{PS}_{M})\), where NF and NM are the total numbers of female and male offspring, respectively, BSF and BSM are the body sizes of female and male offspring, respectively, and PSF and PSM are the survival probabilities of female and male offspring under starvation, respectively, estimated from the Kaplan-Meier survival analysis. The dependent variables in the analyses addressing H3 were the number of females, males, or the total number of offspring. We excluded fleas that did not have any offspring from these analyses. Since any of the interactions between the strain-related variables and age were not significant (Table 2), we repeated the analyses without the age effects. For post-hoc pairwise comparisons, we used the Fisher’s Least Significant Difference test.
Analyses were conducted in R/RStudio44,50, using the packages stats44, agricolae51, and rstatix52.
Results
H1: Two Wolbachia strains coinfect individual S. cleopatrae fleas.
The NGMLST resulted in a total of 2,069,639 reads, distributed similarly among genes (coxA 3399.5 ± 3388.8, ftsZ 3607.3 ± 3482.5, hcpA 2316.7 ± 2082.6, gatB 3605.5 ± 3065.3, fbpA 1614.2 ± 1290.8). In 111 of the 128 tested female fleas, we found two variants for the fbpA, ftsZ, gatB genes and one variant for the coxA and hcpA genes. However, repeating the assays with more generic redesigned primers for the two latter genes resulted in the same variants as for the fbpA, ftsZ, gatB genes.
In the 17 (six from the laboratory and 11 from the field) remaining flea samples, only one variant was found in all the target genes. These 17 samples had identical sequences within genes. According to all phylogenetic trees, this sequence combination was clustered with Wolbachia supergroup B (Fig. 1 and Supplementary Figures S2–S6).
The WGS analysis of the pool of these DNA extracts resulted in a total of 13,632,934 reads and had only an overall 91.5%, identity to Wolbachia wCfeJ. In contrast, it was better matched to the reference Wolbachia sequences clustaring with supergroup B, including Wolbachia strains wDci (overall identity 94.46%), wLst (overall identity 94.55%), and wCon (overall identity 94.9%). We designated this strain as wSc1.
To characterize the second potential coinfecting strain, designated as wSc2, we conducted variant-specific qPCRs on an additional 330 males and female Wolbachia-positive S. cleopatrae fleas from our laboratory colony (see section “qPCRs”). As a result, we detected 47 flea extracts that were only amplified with the wSc2-specific primers. According to all phylogenetic trees, this sequence combination is clustered with Wolbachia wCfeJ (Fig. 1 and Supplementary Figures S2–S6). While the association of this Wolbachia strain with a specific supergroup is not clear, as the different phylogenetic trees reflect slightly different associations with supergroups C, D, and F (Supplementary Figure S2–S6), based on the concatenated tree, both wCfeJ and wSc2 are associated, to some degree, with supergroups C, D, and F, subtending supergroup C (Fig. 1 and Supplementary Figures S2–S6). The WGS analysis of the pool of fleas infected with wSc2 resulted in a total of 6,224,736 reads and had an overall 99.59% identity to wCfeJ. In contrast, it had a lower match to the other reference Wolbachia sequences, namely wDci (overall identity 91.15%), and wTco (overall identity 91.38%). Importantly, the WGS of wSc2 had only a 91.5% identity to the WGS of wSc1.
Unfortunately, because of inadequate Wolbachia coverage in the whole genome sequencing data, which averaged less than 1% of reads, it was not feasible to assemble the Wolbachia genomes into FASTA files and fully complete the wSc1 and wSc2 genomes. Consequently, it was not possible to understand specific gene sequences and functions and utilize them for a more accurate phylogenetic tree.
DNA extracts from a total of 74 Wolbachia-positive female fleas from a previous study that manipulated Wolbachia infection6 were analyzed, both with the Wolbachia-general qPCR and the strain-specific (wSc1 and wSc2) primers. The overall Wolbachia load, as estimated by the sum of the loads calculated by the two variant-specific assays, explained 97% of the variation in the Wolbachia load, as estimated by the general assay (Fig. 2). Moreover, the relationship between the two methods of load estimation was highly significant (F1,72 = 2420, p = 2.2 × 10−16), and the slope was not significantly different from 1 (0.96, 95% CI [0.97, 1.05]), supporting the hypothesis that only Wolbachia wSc1 and wSc2 circulate in the S. cleopatrae flea populations.
Linear correlation between the sum of bacterial load quantified by Wolbachia strain-specific quantitative real-time PCRs (qPCRs) and the bacterial load, as quantified by the Wolbachia-general qPCR assay in 74 extracts of Synosternus cleopatrae fleas. Axes show the number of cells per 5 µl. wSc1, Wolbachia strain wSc1; wSc2, Wolbachia strain wSc2.
H2: Evidence for a fitness advantage provided by Wolbachia wSc1 but not Wolbachia wSc2.
As there was no significant interaction between maternal age with any of the strain-related variables, we repeated the analyses of the integrative index of reproductive success for the entire pool of mother fleas (Table 2). There was no effect of any of the strain loads on the reproductive success index (Table 2). However, after we divided the fleas into five groups, specified based on their naturally low/high loads of wSc1 and wSc2 composition, we found that the reproductive success index was significantly greater in fleas that had high wSc1−low wSc2 than in fleas that had high loads of both strains, low loads of both strains, or had no Wolbachia (Table 2; Fig. 3).
Means ± standard errors (SEs) of the integrated index of reproductive success of female fleas as a function of their Wolbachia strain composition (wSc1 and wSc2). We obtained the strain composition by assigning ‘high’ or ‘low’ status in reference to the median load of each strain. As a baseline comparison, we included in the figure and related analyses the values of the Wolbachia-free female fleas. The letters indicate the Fisher’s Least Significant Difference test results. wSc1, Wolbachia strain wSc1; wSc2, Wolbachia strain wSc2.
H3: Evidence for male killing induced by Wolbachia wSc2 but not by Wolbachia wSc1.
As there was no significant interaction between maternal age with any of the strain-related variables, we repeated the analyses of the number of females, males, and the total number of offspring for the entire pool of mother fleas (Table 2). The number of male offspring was negatively correlated with wSc2 loads (Table 2; Fig. 4A), whereas the number of female offspring was not significantly correlated with any of the strain loads (Table 2; Fig. 4B). Wolbachia wSc2 loads had a similar negative effect on the overall number of offspring (Table 2; Fig. 4C), supporting the mechanism of male killing rather than feminization. Both the strain ratio and the dichotomic wSc1− wSc2 composition were not significantly correlated with the number of male, female, or total offspring (Table 2).
Associations between Wolbachia wSc2 load (number of cells per 5 µl) and the means ± standard errors (SEs) of the numbers of male (A), female (B), and total offspring (C). For visualization purposes, we binned bacterial loads into categories of 1000 cells. wSc2, Wolbachia strain 2.
Discussion
In nature, endosymbiont coinfections are widespread. Yet, they are rarely studied in the context of the endosymbionts that exhibit continuously high incidence with unclear drivers explaining this high prevalence53,54. Our results provide initial evidence that contrasting effects of Wolbachia strains may account for their continuously high incidence in S. cleopatrae flea populations, and similar dynamics may occur in other endosymbiont-arthropod systems. We demonstrated that two Wolbachia strains circulate in natural populations of S. cleopatrae fleas and frequently coinfect them (see also, Supplementary Materials, Fig. S7). We further showed correlative evidence that wSc1 may provide a direct fitness advantage while wSc2 may cause male killing, and hence be costly for fleas. Below we discuss the results supporting the three non-mutually exclusive hypotheses in a broader context, highlight their main insights, and propose future directions to explain the continuously high endosymbiont incidence with unclear drivers in arthropods.
H1: Two Wolbachia strains coinfect S. cleopatrae flea populations.
We characterized two strains of Wolbachia in S. cleopatrae fleas. Given our successful elimination of both Wolbachia strains from fleas6, it is evident that the observed Wolbachia sequences represent genuine strains rather than integrated segments within the host genome. The first strain, designated as wSc1, is clustered with Wolbachia supergroup B. This supergroup consists of Wolbachia bacteria found in a wide range of insects55,56,57. Considering that Wolbachia bacteria can cross host species barriers55, we speculate that at some point during their evolutionary history, wSc1 was acquired by fleas from other insect groups via horizontal transmission.
The second strain, designated as wSc2, is nearly identical to the wCfeJ strain found in C. felis. Based on the MLST analysis, both strains (wSc2 and wCfeJ) are related to supergroups C, D, and F, subtending supergroup C (see also22). The clustering of flea-derived Wolbachia with filarial nematode supergroups58,59 may suggest that these fleas could act as vectors for Wolbachia-infected nematodes (e.g60). However, a prior study of helminth eggs in feces from 78 Gerbillus rodents in the same area found only two non-filarial nematode species—Mastophorus muris and Syphacia sp.—at low prevalence. To further test whether this flea species carries nematodes, we mapped its DNA reads against filarial nematode mitochondrial genomes and found no evidence for the nematode DNA (Supplementary Text S4).
Alternatively, it is possible that in the past, wSc2 was horizontally transmitted from filarial nematodes to fleas, as previously suggested in other systems60,61. To further test the speculations about wSc1 and wSc2’s origin in fleas, additional phylogenomic analyses should be performed, including sequences of closely related bacteria in other host taxa, especially fleas.
H2: Evidence for a fitness advantage provided by Wolbachia wSc1 but not Wolbachia wSc2.
Previously, we did not find significant differences between the current reproductive success of Wolbachia-positive and Wolbachia-free (cured by antibiotic) female fleas6. However, the distinction between the two Wolbachia strains uncovered here revealed that fleas with naturally high wSc1 loads and low wSc2 loads had significantly higher reproductive success than those with high levels of both strains, low levels of both, or no Wolbachia infection. This result suggests that the potential direct advantage posed by wSc1 may be masked by the possibly negative effects of wSc2 on their flea hosts.
According to our phylogenetic analysis, wSc1 clusters with Wolbachia supergroup B, which predominantly consists of reproductive manipulators13. Nevertheless, mutualism and reproductive parasitism by the same strain are not mutually exclusive. Gavotte et al.62 have shown that Wolbachia in Aedes albopictus mosquitoes shift from parasitism to mutualism in response to changes in larval density. Thus, it is possible that wSc1 originally manipulated reproduction in fleas, but since it provided some indirect benefits to female fleas, their relationship has gradually become mutualistic. The nature of relationships between wSc1 and their S. cleopatrae fleas should be further confirmed by (i) comparing the fitness of fleas that are infected solely by wSc1 with those that are Wolbachia-free, and (ii) testing whether the Wolbachia cells previously detected in the flea Malpighian tubules6 are of the wSc1 strain. The latter could indicate that wSc1 plays a nutritional role in fleas.
H3: Evidence for male killing induced by Wolbachia wSc2 but not Wolbachia wSc1.
Female fleas harboring high loads of wSc2 had significantly fewer male and total offspring than females with low loads of wSc2. Interestingly, there was no evidence for sex ratio distortion in Wolbachia-positive versus Wolbachia-free colonies in the previous infection manipulation study by Flatau et al.6. This could indicate that male killing is density-dependent, which was mitigated in the previous study, as the Wolbachia-positive colony consists of females with varying wSc2 loads. The reduction in both male and total offspring observed here suggests the occurrence of male killing rather than feminization induced by Wolbachia. This is because the latter mechanism is not expected to affect the overall offspring number. As in other host species, Wolbachia in S. cleopatrae is transmitted maternally6, suggesting that males are evolutionary dead ends for these endosymbiotic bacteria. Therefore, in theory, male killing should enhance Wolbachia transmission if the fitness of female fleas increases because of the death of their male siblings14. In S. cleopatrae fleas, a reduction in male offspring might reduce the intraspecific competition that is known to occur between larvae63,64, adults65, or both.
Contrary to wSc2, wSc1 loads were not associated with male, female, or overall offspring numbers. Coinfecting Wolbachia strains may affect their host in similar ways but the intensity of their effects may differ66,67. Alternatively, coinfecting Wolbachia strains can play distinct roles in their host biology. For example, it was shown experimentally that the Asobara tabida wasp harbors three strains of Wolbachia, two of which induce CI, while the third one is an obligate endosymbiont necessary to achieve oogenesis68,69. Likewise, genomic analyses suggest that in C. felis fleas, wCfeT is a nutritional endosymbiont, while the coinfecting wCfeJ induces CI22. The nature of the relationships between wSc2 and their S. cleopatrae fleas should be further tested by (i) comparing the offspring number and sex ratio of fleas that are experimentally infected solely with wSc2 with those that are Wolbachia-free and (ii) quantifying larvae and adult intraspecific competition in wSc2-infected offspring that are either supplemented artificially with male fleas or not.
Broad implications of Wolbachia coinfection in flea populations.
This study provides initial indications that the intricate tripartite relationship between two Wolbachia strains and their flea host may be the key toward a better understanding of both the continuously high incidence of Wolbachia in female but not in male fleas, and the high strain coinfection rates observed in natural flea populations (Supplementary Materials, Fig. S7). In particular, we collected evidence suggesting that Wolbachia wSc1 may be responsible for the continuously high incidence observed in female fleas, whereas male killing by Wolbachia wSc2 may explain the low abundance of Wolbachia in the male fleas. Such sex-dependent interactions between a host and its coinfecting endosymbionts may promote the coexistence of the two strains.
The high percentage of flea eggs infected by both Wolbachia strains (Supplementary Fig. S8), along with the stable presence of coinfection observed in both laboratory and natural settings (Supplementary Fig. S7), suggests that these strains are vertically cotransmitted. This adds to the growing evidence of widespread coinfection among vertically transmitted symbionts23,70,71,72. While such transmission typically reduces endosymbiont diversity due to bottlenecks, theoretical models and empirical studies suggest that CI73, functional complementarity74, and interactions between beneficial symbionts and reproductive manipulators75,76,77 can stabilize coinfections over time23. Our findings align with this latter scenario, suggesting a stable long-term coexistence—likely facilitated by cooperation78,79—between the Wolbachia strains wSc1 and wSc2 in natural flea populations. Future studies should explore the nature of associations between these two strains in their flea host and whether they are context-dependent80.
Taken together, our findings contribute to a broader understanding of how endosymbiont diversity is preserved over time despite potential host-symbiont conflicts and changing ecological conditions and highlight the need to pay attention to multiple mechanisms and players under various conditions.
Data availability
All data supporting the findings of this study are available within the paper and its Supplementary Information. MLST raw sequencing data generated in this study have been deposited in the NCBI Sequence Read Archive (SRA) and are available under BioProject (accession number: PRJNA1281168) at https://www.ncbi.nlm.nih.gov/bioproject/1281168.
References
Moran, N. A. Symbiosis Curr. Biology, 16(20): R866–R871. (2006).
Lo, N. et al. Widespread distribution and high prevalence of an alpha-proteobacterial symbiont in the tick Ixodes ricinus. Environ. Microbiol. 8 (7), 1280–1287 (2006).
Douglas, A. E. Nutritional interactions in insect-microbial symbioses: aphids and their symbiotic bacteria buchnera. Ann. Rev. Entomol. 43, 17–37 (1998).
Hosokawa, T. et al. Infection prevalence of Sodalis symbionts among Stinkbugs. Zoological Lett. 1, 1–7 (2015).
Richardson, K. M. et al. A Wolbachia infection from Drosophila that causes cytoplasmic incompatibility despite low prevalence and densities in males. Heredity 122 (4), 428–440 (2019).
Flatau, R., Segoli, M. & Hawlena, H. Wolbachia endosymbionts of fleas occur in all females but rarely in males and do not show evidence of obligatory relationships, fitness effects, or sex-distorting manipulations. Front. Microbiol., 12, ARTN 649248 (2021).
Ferrari, J. & Vavre, F. Bacterial symbionts in insects or the story of communities affecting communities. Philosophical Trans. Royal Soc. B-Biological Sci. 366 (1569), 1389–1400 (2011).
Feldhaar, H. Bacterial symbionts as mediators of ecologically important traits of insect hosts. Ecol. Entomol. 36 (5), 533–543 (2011).
Su, Q., Zhou, X. & Zhang, Y. Symbiont-mediated functions in insect hosts. Commun. Integr. Biol., 6, ARTN e23804 (2013).
Cao, L. J., Jiang, W. B. & Hoffmann, A. A. Life history effects linked to an advantage for wAu Wolbachia in Drosophila. Insects, 10(5), ARTN 126 (2019).
Brownlie, J. C. et al. Evidence for metabolic provisioning by a common invertebrate endosymbiont, Wolbachia pipientis, during periods of nutritional stress. PLoS Pathog. 5 (4), ARTNe1000368 (2009).
Li, C. F. et al. Co-infection with Wolbachia and cardinium May promote the synthesis of fat and free amino acids in a small spider, Hylyphantes graminicola. J. Invertebr. Pathol. 169, ARTN107307 (2020).
Werren, J. H., Baldo, L. & Clark, M. E. Wolbachia: master manipulators of invertebrate biology. Nat. Rev. Microbiol. 6 (10), 741–751 (2008).
Hurst, G. D. D. & Jiggins, F. M. Male-killing bacteria in insects: mechanisms, incidence, and implications. Emerg. Infect. Dis. 6 (4), 329–336 (2000).
Bouchon, D., Cordaux, R. & Grève, P. Feminizing Wolbachia and the evolution of sex determination in isopods. Insect Symbiosis. 3, 273–294 (2008).
Engelstädter, J. & Hurst, G. D. D. The ecology and evolution of microbes that manipulate host reproduction. Annu. Rev. Ecol. Evol. Syst. 40, 127–149 (2009).
Wedincamp, J. Jr. & Foil, L. D. Vertical transmission of Rickettsia felis in the Cat flea (Ctenocephalides felis Bouche). J. Vector Ecol. 27 (1), 96–101 (2002).
Charlat, S. et al. Prevalence and penetrance variation of male-killing Wolbachia across Indo-Pacific populations of the butterfly Hypolimnas Bolina. Mol. Ecol. 14 (11), 3525–3530 (2005).
White, J. A. et al. Endosymbiont costs and benefits in a parasitoid infected with both Wolbachia and cardinium. Heredity 106 (4), 585–591 (2011).
Hamm, C. A. et al. Wolbachia do not live by reproductive manipulation alone: infection polymorphism in Drosophila Suzukii and D. subpulchrella. Mol. Ecol. 23 (19), 4871–4885 (2014).
Flatau, R. et al. Wolbachia’s role in mediating its flea’s reproductive success differs according to flea origin. FEMS Microbiol. Ecol. 94 (10), ARTNfiy157 (2018).
Driscoll, T. P. et al. Evolution of Wolbachia mutualism and reproductive parasitism: insight from two novel strains that co-infect Cat fleas. Peerj 8, ARTNe10646 (2020).
Vautrin, E. & Vavre, F. Interactions between vertically transmitted symbionts: Cooperation or conflict? Trends Microbiol. 17 (3), 95–99 (2009).
Kikuchi, Y. & Fukatsu, T. Diversity of Wolbachia endosymbionts in heteropteran Bugs. Appl. Environ. Microbiol. 69 (10), 6082–6090 (2003).
Flórez, L. V. & Kaltenpoth, M. Symbiont dynamics and strain diversity in the defensive mutualism between Lagria beetles and Burkholderia. Environ. Microbiol. 19 (9), 3674–3688 (2017).
McLean, A. H. C. et al. Consequences of symbiont co-infections for insect host phenotypes. J. Anim. Ecol. 87 (2), 478–488 (2018).
Wei, J. et al. High prevalence of Wolbachia infection does not explain unidirectional cytoplasmic incompatibility of Altica flea beetles. Ecol. Entomol. 45 (1), 67–78 (2020).
Halter, T. et al. One to host them all: genomics of the diverse bacterial endosymbionts of the spider Oedothorax gibbosus. Microb. Genomics. 9 (2), ARTNmgen000943 (2023).
Richardson, K. M. et al. A male-killing Wolbachia endosymbiont is concealed by another endosymbiont and a nuclear suppressor. PLoS Biol. 21 (3), ARTNe3001879 (2023).
Ant, T. H. & Sinkins, S. P. A Wolbachia triple-strain infection generates self-incompatibility in Aedes albopictus and transmission instability in Aedes aegypti. Parasites Vectors, 11(1), ARTN 295 (2018).
Arai, H. et al. Multiple infection and reproductive manipulations of Wolbachia in Homona magnanima in (Lepidoptera: Tortricidae). Microb. Ecol. 77 (1), 257–266 (2019).
Dedeine, F. et al. Intra-individual coexistence of a Wolbachia strain required for host oogenesis with two strains inducing cytoplasmic incompatibility in the Wasp Asobara tabida. Evolution 58 (10), 2167–2174 (2004).
Jones, M. W. et al. Infection dynamics of cotransmitted reproductive symbionts are mediated by sex, tissue, and development. Appl. Environ. Microbiol. 88 (13), ARTNe00529–ARTNe00522 (2022).
Hiroki, M. et al. Multiple infection with Wolbachia inducing different reproductive manipulations in the butterfly Eurema hecabe. Proc. R. Soc. Lond. B Biol. Sci. 271 (1549), 1751–1755 (2004).
Łukasik, P. et al. Unrelated facultative endosymbionts protect aphids against a fungal pathogen. Ecol. Lett. 16 (2), 214–218 (2013).
Asselin, A. K. et al. Contrasting patterns of virus protection and functional incompatibility genes in two conspecific Wolbachia strains from Drosophila Pandora. Appl. Environ. Microbiol. 85 (5), ARTNe02290–ARTNe02218 (2019).
Liang, X. et al. Wolbachia Inter-Strain competition and Inhibition of expression of cytoplasmic incompatibility in mosquito. Front. Microbiol. 11, ARTN1638 (2020).
Narita, S., Nomura, M. & Kageyama, D. Naturally occurring single and double infection with Wolbachia strains in the butterfly Eurema hecabe: transmission efficiencies and population density dynamics of each Wolbachia strain. FEMS Microbiol. Ecol. 61 (2), 235–245 (2007).
Watanabe, M. et al. Superinfection of cytoplasmic incompatibility-inducing Wolbachia is not additive in Orius strigicollis (Hemiptera: Anthocoridae). Heredity 106 (4), 642–648 (2011).
Curry, M. M. et al. Multiple endosymbiont infections and reproductive manipulations in a linyphiid spider population. Heredity 115 (2), 146–152 (2015).
Werren, J. H., Skinner, S. W. & Huger, A. M. Male-killing bacteria in a parasitic Wasp. Science 231 (4741), 990–992 (1986).
Messika, I. et al. From endosymbionts to host communities: factors determining the reproductive success of arthropod vectors. Oecologia 184 (4), 859–871 (2017).
Baldo, L. et al. Multilocus sequence typing system for the endosymbiont Wolbachia pipientis. Appl. Environ. Microbiol. 72 (11), 7098–7110 (2006).
Core Team, R. R., R: A Language and Environment for Statistical Computing (R Foundation for Statistical Computing, 2022).
Callahan, B. J. et al. DADA2: High-resolution sample inference from illumina amplicon data. Nat. Methods. 13 (7), 581–583 (2016).
Sharma, A. K. & Som, A. Assigning new supergroups V and W to the Wolbachia diversity. Bioinformation 19 (3), 336–340 (2023).
Guindon, S. et al. New algorithms and methods to estimate maximum-likelihood phylogenies: assessing the performance of PhyML 3.0. Syst. Biol. 59 (3), 307–321 (2010).
Baneshi, M. R. & Talei, A. R. Dichotomisation of continuous data: review of methods, advantages, and disadvantages. Iran. J. Cancer Prev. 4 (1), 26–32 (2011).
Williamson, S. F. et al. Subgroup analyses in randomized controlled trials frequently categorized continuous subgroup information. J. Clin. Epidemiol. 150, 72–79 (2022).
Team, R. RStudio: Integrated development environment for R. Posit Software, PBC: Boston, MA. (2022).
de Mendiburu, F. Agricolae Tutorial (Version 1.3-5) (La Molina, Peru, 2021).
Kassambara, A. Pipe-friendly Framework for Basic Statistical Tests [R Package Rstatix Version 0.7. 2] (Vienna, Austria, 2023). Comprehensive R.
Kanthong, N. et al. Persistent, triple-virus co-infections in mosquito cells. BMC Microbiol. 10, ARTN14 (2010).
Monack, D. M., Mueller, A. & Falkow, S. Persistent bacterial infections: the interface of the pathogen and the host immune system. Nat. Rev. Microbiol. 2 (9), 747–765 (2004).
Raychoudhury, R. et al. Modes of acquisition of wolbachia: horizontal transfer, hybrid introgression, and codivergence in the Nasonia species complex. Evolution 63 (1), 165–183 (2009).
Shaikevich, E. et al. Wolbachia symbionts in mosquitoes: Intra- and intersupergroup recombinations, horizontal transmission and evolution. Mol. Phylogenet. Evol. 134, 24–34 (2019).
Vancaester, E. & Blaxter, M. Phylogenomic analysis of Wolbachia genomes from the Darwin tree of life biodiversity genomics project. PLoS Biol., 21(1), ARTN e3001972 (2023).
Ferri, E. et al. New insights into the evolution of Wolbachia infections in filarial nematodes inferred from a large range of screened species. Plos One. 6 (6), p–ARTNe20843 (2011).
Lefoulon, E. et al. Pseudoscorpion Wolbachia symbionts: diversity and evidence for a new supergroup S. BMC Microbiol. 20 (1), ARTN188 (2020).
Traversa, D. Fleas infesting pets in the era of emerging extra-intestinal nematodes. Parasites Vectors. 6, ARTN59 (2013).
Vavre, F. et al. Phylogenetic evidence for horizontal transmission of Wolbachia in host-parasitoid associations. Mol. Biol. Evol. 16 (12), 1711–1723 (1999).
Gavotte, L. et al. Costs and benefits of Wolbachia infection in immature Aedes albopictus depend upon sex and competition level. J. Invertebr. Pathol. 105 (3), 341–346 (2010).
Lawrence, W. & Foil, L. D. The effects of diet upon pupal development and cocoon formation by the Cat flea (Siphonaptera: Pulicidae). J. Vector Ecol. 27 (1), 39–43 (2002).
Khokhlova, I. S. et al. Adaptation to a novel host and performance trade-off in host-generalist and host-specific insect ectoparasites. Insect Sci. 29 (2), 567–580 (2022).
Hawlena, H. et al. Host defence versus intraspecific competition in the regulation of infrapopulations of the flea Xenopsylla conformis on its rodent host Meriones crassus. Int. J. Parasitol. 37 (8–9), 919–925 (2007).
Kondo, N. et al. Prevailing triple infection with Wolbachia in Callosobruchus chinensis (Coleoptera: Bruchidae). Mol. Ecol. 11 (2), 167–180 (2002).
Kondo, N. et al. Genome fragment of Wolbachia endosymbiont transferred to X chromosome of host insect. Proc. Natl. Acad. Sci. U.S.A. 99 (22), 14280–14285 (2002).
Mouton, L. et al. Virulence, multiple infections and regulation of symbiotic population in the Wolbachia-Asobara Tabida symbiosis. Genetics 168 (1), 181–189 (2004).
Dedeine, F. et al. Removing symbiotic Wolbachia bacteria specifically inhibits oogenesis in a parasitic Wasp. Proc. Natl. Acad. Sci. U.S.A. 98 (11), 6247–6252 (2001).
Zhao, D. et al. Interactions between facultative symbionts Hamiltonella and cardinium in Bemisia tabaci (Hemiptera: Aleyrodoidea): cooperation or conflict? J. Econ. Entomol. 111 (6), 2660–2666 (2018).
Zhu, Y. X. et al. Spider mites singly infected with either Wolbachia or Spiroplasma have reduced thermal tolerance. Front. Microbiol., 12, ARTN 706321 (2021).
Ueda, M. et al. Distinct effects of three Wolbachia strains on fitness and immune traits in Homona magnanima. Heredity 130 (1), 22–29 (2023).
Mouton, L. et al. Strain-specific regulation of intracellular Wolbachia density in multiply infected insects. Mol. Ecol. 12 (12), 3459–3465 (2003).
Buysse, M. et al. A dual endosymbiosis supports nutritional adaptation to hematophagy in the invasive tick Hyalomma marginatum. Elife 10, ARTNe72747 (2021).
Tsuchida, T. et al. Diversity and geographic distribution of secondary endosymbiotic bacteria in natural populations of the pea aphid, Acyrthosiphon pisum. Mol. Ecol. 11 (10), 2123–2135 (2002).
Zchori-Fein, E. & Brown, J. K. Diversity of prokaryotes associated with Bemisia tabaci (Gennadius) (Hemiptera: Aleyrodidae). Ann. Entomol. Soc. Am. 95 (6), 711–718 (2002).
Gomez-Valero, L. et al. Coexistence of Wolbachia with Buchnera aphidicola and a secondary symbiont in the aphid Cinara cedri. J. Bacteriol. 186 (19), 6626–6633 (2004).
Vautrin, E. et al. Do vertically transmitted symbionts co-existing in a single host compete or cooperate? A modelling approach. J. Evol. Biol. 21 (1), 145–161 (2008).
Zélé, F. et al. Ecology and evolution of facilitation among symbionts. Nat. Commun., 9, ARTN 4869 (2018).
Rock, D. I. et al. Context-dependent vertical transmission shapes strong endosymbiont community structure in the pea aphid, Acyrthosiphon pisum. Mol. Ecol. 27 (8), 2039–2056 (2018).
Acknowledgements
This study was supported by an Israel Science Foundation grant (award 1391/15 to HH) and an Ecology of Infectious Diseases grant (award DEB-1813069 from the National Science Foundation) under the auspices of the US-Israel Binational Science Foundation to HH and JEB. AK was supported by the Excellence Fellowship Program for International Postdoctoral Researchers from The Israel Academy of Sciences and Humanities and by a fellowship from the Kreitman School of Advanced Graduate Studies at Ben-Gurion University of the Negev. We thank Daniel E. Deatherage for assistance with next-generation sequencing. Arthropod and nematode icons used in Figure 1 and Supplementary Figures S2–S6 were downloaded from open access databases: PNGegg and Silhouette AC or created based on the graphics from Driscoll et al. (Figure 4;22).
Author information
Authors and Affiliations
Contributions
Project Administration: HH; Resources: HH; RF, MS, and HH conceived the study. RF and AK collected the data. RF performed molecular analyses. JB analyzed WGS data. RF, AK, and HH conducted statistical analyses. RF, AK, MS, and HH interpreted the results. AK visualized the results. AK and HH wrote the first version of the manuscript. RF, MS, JB, and HH revised, edited, and approved the manuscript.
Corresponding author
Ethics declarations
Competing interests
The authors declare no competing interests.
Additional information
Publisher’s note
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Electronic supplementary material
Below is the link to the electronic supplementary material.
Rights and permissions
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.
About this article
Cite this article
Flatau, R., Krawczyk, A.I., Segoli, M. et al. Continuously high Wolbachia incidence in flea populations may result from dual-strain infections with divergent effects. Sci Rep 15, 21720 (2025). https://doi.org/10.1038/s41598-025-09403-2
Received:
Accepted:
Published:
DOI: https://doi.org/10.1038/s41598-025-09403-2
