The Wolbachia mobilome in Culex pipiens includes a putative plasmid

Wolbachia is a genus of obligate intracellular bacteria found in nematodes and arthropods worldwide, including insect vectors that transmit dengue, West Nile, and Zika viruses. Wolbachia’s unique ability to alter host reproductive behavior through its temperate bacteriophage WO has enabled the development of new vector control strategies. However, our understanding of Wolbachia’s mobilome beyond its bacteriophages is incomplete. Here, we reconstruct near-complete Wolbachia genomes from individual ovary metagenomes of four wild Culex pipiens mosquitoes captured in France. In addition to viral genes missing from the Wolbachia reference genome, we identify a putative plasmid (pWCP), consisting of a 9.23-kbp circular element with 14 genes. We validate its presence in additional Culex pipiens mosquitoes using PCR, long-read sequencing, and screening of existing metagenomes. The discovery of this previously unrecognized extrachromosomal element opens additional possibilities for genetic manipulation of Wolbachia.

M osquitoes are major vectors of disease-causing pathogens worldwide including viruses such as dengue, West Nile, Chikungunya, Zika, and yellow fever [1][2][3] . In the absence of effective vaccines and the off-target effects of insecticides used to control mosquitoes, novel vector biocontrol efforts are the focus of intense study 4 . Over the past decade, the widely distributed endosymbiotic alphaproteobacteria Wolbachia has gained attention as a promising mosquito-control strategy. The basic reasons for its heightened attention is that Wolbachia can cause reproductive parasitism whereby the intracellular bacteria in the reproductive tissues can alter sexual reproduction to enhance its maternal spread through host populations at the expense of host fitness [5][6][7] . Second and notably, native Wolbachia reduce RNA virus replication in the fruit fly Drosophila 8-10 , and Wolbachia transinfections into mosquito species (i.e., Aedes aegypti, Ae. albopictus, Ae. polynesienses, and Anopheles stephensi) similarly result in mosquito lines refractory to various types of pathogenic RNA viruses [11][12][13][14] . The combination of a selfish drive system and pathogen blocking by Wolbachia has led to successful pilot trials for suppression of mosquito population size and replacement of infected mosquito populations so they can no longer transmit pathogens [15][16][17] .
Wolbachia are transovarially transmitted from the mother to offspring 6,[18][19][20] . In some arthropods, including the naturally infected vector species complex Culex pipiens, Wolbachia 'modify' sperm in testes, leading to embryonic lethality if the infected male mates with either an uninfected female or a female harbouring an incompatible Wolbachia strain. When both male and female are infected with the same Wolbachia, the modification is 'rescued', and compatibility is restored 21 . This reproductive alteration, termed 'cytoplasmic incompatibility' (CI), comprises the most common form of Wolbachia-induced parasitism and, as studied in C. pipiens, can lead to highly diverse unidirectional and bidirectional incompatibility phenotypes 22 . CI is currently used in vector control studies for population replacement by Wolbachiainfected strains that block arbovirus transmission 17 or population suppression that reduces the number of mosquito vectors 23 . Notably, the genes responsible for CI [24][25][26] occur in the eukaryotic association module of Wolbachia's temperate phage WO [25][26][27] , which laterally transfers between Wolbachia coinfections and evolves rapidly [28][29][30][31][32] . Overall, these findings emphasize the importance of further investigating mobile genetic elements in Wolbachia.
While recent studies shed light on the role of phage WO in Wolbachia genome evolution and CI, other extrachromosomal elements, such as plasmids, have not been detected in the symbiont. Notably, over half of the species in the closely related Rickettsia genus have plasmids 33 that play roles in DNA replication, partitioning, mobilization, and conjugation 34,35 and offer a potential tool for genetic manipulation of diverse members of Rickettsia 36,37 . Similar genetic manipulation strategies for Wolbachia are conceivable 27 ; however, previous efforts to search for such extrachromosomal mobile genetic elements have not been successful 38,39 . The lack of isolates limit direct insights into Wolbachia genomics, and most metagenomic approaches thus far rely on pooled individuals grown in the laboratory environment due to low infection densities 40,41 . Since these limitations can conceal naturally occurring genomic diversity among Wolbachia populations, highly resolved analyses of individual mosquitoes may reveal additional insights into the Wolbachia 'mobilome', the pool of all mobile genetic elements associated with Wolbachia populations.
Here we sequence ovary samples from four wild-caught C. pipiens individuals captured in Southern France from a single trap. Using genome-resolved metagenomic and pangenomic analysis strategies, we were able to reconstruct and compare nearcomplete Wolbachia genomes from each individual. Besides a diverse set of virus-associated genes that were missing or absent in the reference Wolbachia genome wPip Pel, our data reveal the first lines of evidence for an extrachromosomal circular element with genetic and functional hallmarks of a plasmid that we tentatively name pWCP.

Results
A Wolbachia metagenome-assembled genome (MAG) is recovered in each sample. Shotgun sequencing of DNA recovered from ovary samples of four C. pipiens individuals (O03, O07, O11, O12) resulted in 65-78 million paired-end sequences after quality filtering. Metagenomic assembly of each sample individually yielded 147K-183K contiguous DNA segments (contigs) >1 kbp, which recruited 48.1-72.9% of the raw sequencing reads. The relatively high fraction of unmapped reads were likely due to challenges associated with the assembly of environmental metagenomes 42,43 , especially in the presence of eukaryotic host genomic DNA. Supplementary Table 1 reports statistics for the raw number of reads and assembly results for each sample. We employed a metagenomic binning strategy that uses sequence composition signatures and differential coverage statistics of contigs across samples. For each ovary sample, we were able to reconstruct a highly complete single bacterial MAG that resolved to Wolbachia (Table 1).
Wolbachia MAGs include highly covered non-phage contigs. The relatively low number of single-nucleotide variants (0.01-0.05%, Supplementary pWCP: a Wolbachia-associated putative plasmid. Based on homology between the ends of these five contigs and the flanking regions of IS110 TE, we predicted that the missing region was the latter element. We artificially circularized Contig O11_A (8037 bp) and inserted an IS110 TE (1386 bp) based on the overlapping 5' and 3' ends ( Fig. 1). Metagenomic read recruitment onto this artificially circularized contig, which we tentatively name 'pWCP' (for plasmid of Wolbachia endosymbiont in C. pipiens), showed consistent coverage over its entire length except a clear two-fold coverage increase in a region that matched to the IS110 TE in all four C. pipiens individuals in our study (Supplementary Figure 1). These data suggested the IS110 TE are located in the extrachromosomal pWCP while some others could be integrated in the bacterial chromosome. The read recruitment from three C. pipiens egg metagenomes generated in a previous study 32 confirmed the near identical presence of pWCP in all three, including the increase in coverage matching to the IS110 TE (Supplementary Figure 1). To validate the circular and extrachromosomal nature of pWCP independently of short-read recruitment-and assemblybased strategies, we generated long reads from additional C. pipiens complex samples using a MinION sequencer. Since MinION sequencing occurs with no polymerase chain reaction (PCR) or downstream assembly steps, we hypothesized that long reads that match to pWCP should never be (1) flanked by genomic regions matching to the Wolbachia chromosome and (2) longer than the pWCP itself. Our MinION sequencing analysis resulted in 14,808 high-quality sequences that were >5000 nucleotides. While a significant fraction of these reads were eukaryotic contamination and the lambda phage DNA that we used to pad our low-biomass samples (see Methods), a local BLAST search of artificially circularized pWCP sequence against long reads revealed 19 that aligned to pWCP with an e value of <1e−20 (Supplementary Data 2). Thirteen of these 19 reads covered >50% of the length of pWCP and contained no other genomic region. In other words, each of the long reads were equal to or shorter than the length of pWCP, as expected for an extrachromosomal element (Fig. 2a). Moreover, 6 of these 13 reads were exactly equal to the length of pWCP, covering its entirety with non-identical start positions, confirming its circularity (Fig. 2a). The final 6 long reads that covered <50% of the length of pWCP had specific matches to the IS110 TE (Fig. 2b); this result is consistent with the multiple occurrences of the IS110 TE in the Wolbachia genome and explains the increase in coverage in metagenomic read recruitment results (Supplementary Figure 1).
To further validate the extrachromosomal nature of pWCP, we PCR screened genomic DNA from the wild caught individuals as well as from tetracycline (TC)-treated laboratory lines (TC1-TC3; negative controls). TC treatment eliminates Wolbachia from its hosts and is commonly used to generate C. pipiens uninfected laboratory lines 47 . First, using LCO1490 and HCO2198 mitochondrial cytochrome c oxidase subunit I (COI) invertebrate primers 48 , we detected the presence of a~708-bp band in the four C. pipiens individuals and the three Wolbachia-free Culex quinquefasciatus controls treated with TC (TC1-TC3), confirming the presence of arthropod DNA in all samples (Fig. 1b). Next, we verified the presence of Wolbachia DNA in the first four samples by PCR amplifying a 438-bp fragment of the Wolbachia 16S rRNA gene using Wspec-F and Wspec-R primers 49 . No band was observed in the TC samples, confirming the absence of Wolbachia (Fig. 1c). Critically, a~1800-bp fragment amplified with primers 263F and 2127R (Supplementary Table 5) designed from the ends of the artificially circularized contig within DnaB gene (and uniquely matching those sites), confirmed the circular nature of the pWCP and its presence only in Wolbachia-infected C. pipiens samples (Fig. 1d). We also designed primers to Sanger sequence across the circular gap (see Supplementary Table 5 for primer sets EC_1-EC_7), and results confirmed the pWCP sequence obtained with both Illumina and MinION sequencing. Finally, we amplified IS110 TE with primers EC_4F and EC_4R (Supplementary Table 5) in the four C. pipiens samples studied herein while observing no band in the Wolbachia-free samples (Fig. 1e) in order to verify the strict association of IS110 TE with Wolbachia (Fig. 1f). These results were further confirmed using additional rifampicin-and oxytetracycline-treated C. pipiens samples (Supplementary Figure 2). Of note, we verified the presence of pWCP in the originating TC1-TC3 and C. pipiens laboratory stocks prior to antibiotic treatment.
The average nucleotide identity of the four independently assembled pWCP sequences from individual mosquitoes was 99.65-100% (Supplementary Table 6). In addition, we identified a 8315-bp contig in the wPip JHB assembly (ABZA01000008.1 50 ), which also was >99.53% identical to each of the four pWCP sequences (Supplementary Table 6). The IS110 TE was 100% identical across samples and confirmed as clonal. Our additional read recruitment analyses from publicly available metagenomes (SRR5810516, SRR5810517, SRR5810518) 32 also revealed the widespread occurrence of pWCP in C. pipiens individuals from Turkey, Algeria, and Tunisia (Supplementary Note 1, Supplementary Figure 3) at a similar 4.2-7.3-fold higher copy number relative to the Wolbachia genome. We also confirmed that pWCP and the Wolbachia genome display a similar tetranucleotide composition (Supplementary Note 2, Supplementary Figure 4). Overall, these findings suggest that pWCP is maintained over the long term with Wolbachia from these strains.
pWCP encodes an IS110 TE, 14 genes, and an intergenic repeat region. The IS element is homologous to ISWpi12 of the IS110 family 46 based on IS finder platform search. Annotation of genes in pWCP also revealed the presence of a disrupted DnaB-like helicase, two RelBE loci, a ParA-like gene (each with identical AA sequences across samples, Supplementary Note 3), and multiple genes encoding hypothetical proteins (Fig. 1f, Supplementary Data 3). The ParA-like partitioning sequence showed amino acid homology to the chromosome partitioning of plasmid protein (ParA) identified in Ca. Caedibacter acanthamoebae (an endosymbiont of acanthamoebae; e value: 2 × 10 -46 ), the bacterium Odyssella thessalonicensis (e value: 9 × 10 -45 ), and Rickettsia raoultii (e value: 7 × 10 -41 ). The RelBE toxin-antitoxin (TA) locus has been identified in multiple Wolbachia genomes 51 and is often associated with prophage WO regions (e.g., of wVitA, wHa, wMel, wAu, wRi, wSuzi, wFol, wInc), specifically within the tail and/or capsid modules. In other bacteria, this TA system can promote the stability of its encoding mobile element, including plasmids or pathogenicity islands, through post-segregational killing of cells that have lost the antitoxin component of the TA operon 52,53 . Most remaining pWCP genes were hypothetical and unique to either wPip and/or the B-Wolbachia phyletic supergroup 54 (Supplementary Data 3), including GP-09 and GP-11 which showed a very weak homology to a Transcription Factor and a Terminase, respectively (e value > 0 in SCOPe, Pfam, and Clusters of Orthologous Group (COG), highlighting the need for further functional characterization of this newly discovered mobile genetic element. In particular, terminase proteins bind and package DNA into the capsids of phage particles 55 . These data indicate that the extrachromosomal DNA could potentially fall into three categories: a simple plasmid, a mini-chromosome of Wolbachia, or a plasmid-like replicon that hijacks the capsids of phage WO.
Beyond the putative coding regions, alignment of all contigs revealed an intergenic variable number tandem repeat (VNTR) region ( Fig. 3a repeats adjacent to parA that vary in number among individual pWCP sequences. Recent studies observed the same repeat region, identified as pp-hC1A_5, and used it to genotype different strains of Wolbachia 56,57 , yet these have not been studied at an individual level. The authors suggested that a deletion in the intergenic polymorphic region could serve as a recombination or horizontal gene transfer site 56 . Alternatively, we hypothesize that the direct repeats, as present in iteron plasmids, could indicate a potential origin of replication and play a role in copy number control 58 . Our analysis of the pWCP sequence also revealed a  (Fig. 3c). Although the role of these sequences is not clear, the closely related plasmid of Rickettsia monacensis (pRM) harbours four perfect and four imperfect palindromes 34 .
The Wolbachia metapangenome reveals novel viral genes. The assembly and binning of individual mosquitoes from the wild also enabled comparison of the diversity and the gene content of prophage WO regions in our Wolbachia genomes vs. the wPip Pel reference genome. We performed a metapangenomic analysis of the four Wolbachia MAGs and wPip Pel in conjunction with the four metagenomes from individual mosquitoes. By linking genes to their abundances in C. pipiens metagenomes, we aimed at tying genomic and environmental data. To determine gene coverages, we used the Wolbachia MAG O07 as reference for read recruitment since (1) it was the largest MAG in size with most number of genes (Table 1) and (2) all MAGs were over 99.8% identitical (Supplementary Table 3).
The Wolbachia pangenome contained 1166 gene clusters (that is, groups of homologous predicted open reading frames (ORFs) based on amino acid sequence identity), the majority of which were conserved across all five genomes (Fig. 4, Supplementary   Figure 6). wPip Pel and Wolbachia MAG O07 carried the largest number of unique gene clusters (Supplementary Data 4). Genes that were unique to wPip Pel (n = 41) encoded functions including several transposases, bacteriophage capsid protein coding genes, and other phage-related sequences, most of which were associated with known Wolbachia prophages.
Gene clusters unique to MAG O07 (n = 56) included a gene coding for an ankyrin and tetratricopeptide repeat protein previously identified in phage WO from Nasonia vitripennis wasps 27 . Ankyrin and tetratricopeptide repeats mediate a broad range of protein-protein interactions (apoptosis, cell signaling, inflammatory response, etc.) within eukaryotic cells and are commonly associated with effector proteins of certain intracellular pathogens 59,60 . There was also a Retron-type reverse transcriptase and genes coding for Transposases (COG3293 and a Transposase InsO and inactivated derivatives gene). Although most remaining unique O07 gene clusters had no functional annotation, about a third matched to eukaryotic viral genes based on homology searches in the NCBI's non-redundant protein sequence database (Supplementary Data 5).
These data add to previous studies showing that regions of genomic diversity between closely related Wolbachia genomes are often virus associated 50,61-63 . Note that most gene clusters with genes that were unique to MAG O07 did recruit reads from the three other metagenomes (Fig. 4, Supplementary Figure 6), suggesting that, even though they were not characterized in our other MAGs, they did occur in C. pipiens metagenomes. Absence of these genes from our MAGs are most likely due to (1) assembly artifacts that result in fragmented contigs that are too short to be considered for binning or (2) mutations in the gene context that affect the gene caller to identify them properly. Gene clusters that matched to pWCP also occurred only in our Wolbachia MAGs and were missing in wPip Pel (Supplementary Figure 6), which is expected since the wPip Pel reference genome is solely composed of the Wolbachia bacterial chromosome 40 . Overall, the metapangenome sheds light on a substantial amount of viral genetic diversity, revealing almost as many virus-associated genes as the ones that were previously recognized in the reference genome wPip.
Unlike our Wolbachia MAGs, wPip Pel is a high-quality genome assembled into a single scaffold. Even though wPip Pel is not completely closed 40 , it offers a more complete representation of the synteny of genes in comparison to our MAGs. Hence, in addition to ordering gene clusters based on their distribution patterns across all genomes (Supplementary Figure 6), we took advantage of the wPip Pel genome to determine the order of gene clusters in the metapangenome based on the order of wPip genes in the wPip Pel genome (Fig. 4). This 'forced synteny' allowed us to investigate the diversity and abundance of phage genes in the context of the five previously identified prophage WO regions in wPip Pel (Fig. 4). Some gene clusters within prophage regions appeared to be unique to wPip Pel and were not detected in our metagenomes. It is possible that these genes were not recovered from our MAGs due to small contig size (that is, contigs were too short to be considered for binning). However, our MAGs often carried upstream and downstream phage genes in these regions, suggesting that while some phage genes were conserved across all genomes, others differed significantly from their wPip Pel counterparts (Fig. 4, Supplementary Data 4). It is possible that a set of new phage-associated genes only found in MAGs (Fig. 4, Supplementary Data 5) have functional homologues in wPip Pel. Previous studies indeed show WO genes that have distinct nucleotide sequences yet similar predicted functions 30 . However, it is also possible that some genes detected only in our our MAGs at the sequence level may also be encoding unique functions compared to known phage genes; for instance, eukaryotic-like homologues were recently shown as constituents of phage WO 27 .  Metagenomic read recruitment revealed between a 1.5-and 5fold increase in coverage between pWCP and some structural phage genes (e.g., WP0415 and WP0446 Tail genes) in our metagenomes compared to the coverage of the bacterial chromosome in all four C. pipiens individuals (Supplementary Data 6). The forced synteny organization of gene clusters also revealed that a single prophage WO region could include both high-and low-coverage phage genes (Fig. 4). The multi-copy occurrence of pWCP could explain its increased coverage (Supplementary Figure 1), and differential coverage regimes for genes within a single prophage region could be explained by at least two different models. First, increased coverage of some prophage genes could be attributed to lytic activity: the prophage genes displaying lower bacterial-like coverage are not part of the virion, while those with higher coverage correspond to phage genes that are replicated and packaged into phage capsids. This lytic model is consistent with the observation of phage particles in C. pipiens mosquitoes 44,45 , observed lytic activity of Wolbachia phages in Nasonia testes, and the sequencing of WO genomes from purified phage particles 27 . The partial replication and packaging of prophage WO genes could result from either a 'less than headful' mechanism of packaging, as described in model  phages P1 and T4 64,65 , or it could represent active vs. degenerate prophage variants in the wPip Pel chromosome. Second, some genes in prophage WO genomes could be copied throughout the Wolbachia chromosome explaining the increase in coverage due to their multi-copy occurrence. This model is supported by the prophage duplication events observed in the wRi and wSuz genomes 61,66 as well as the presence of TEs within and/or flanking prophage variants 27,40 that could enable genomic rearrangement and duplication. These models may not be mutually exclusive, and the system could involve both duplication of prophage genomes and the induction of phage particles.

Discussion
Shotgun metagenomes from individual C. pipiens ovary samples allowed us to de novo reconstruct Wolbachia genomes from single mosquitoes and compare these MAGs to each other as well as to the reference wPip Pel genome through pangenomic strategies. Our data reveal an extensive diversity of previously undetected Wolbachia phage WO and other viral genes and notably the first indication that Wolbachia harbours a candidate plasmid, shedding new light on the richness of the Wolbachia mobilome. The definition of a plasmid varies; here we adhere to the typical characterization of a plasmid as a small hereditary, extrachromosomal, circular element.
Even though a Wolbachia plasmid has not been reported before, we did find evidence of its occurrence in C. pipiens metagenomes in previous studies. It is likely that the previous efforts overlooked this element due to computational challenges associated with the assembly of metagenomic short reads. Beyond the IS110 TE that occurred both in the plasmid and Wolbachia chromosome, pWCP contained a region of intergenic VNTR that differed across individuals in our study. This suggests that the co-assembly of pooled individuals may yield fragmented contigs. The same VNTR sequences were found in C. pipiens from Australia, Argentina, USA, Italy, Japan, Israel, and Greece 56 , and we confirmed the presence of pWCP in C. pipiens samples from countries of the Mediterranean basin 32 . The high similarity of pWCP sequences across C. pipiens metagenomes in addition to its global prevalence suggest evolutionary constraints and a possible functional role in Wolbachia symbiosis. However, we did not detect the plasmid in other Wolbachia strains through screening available metagenomes from the fly Drosophila melanogaster, the planthopper Laodelphax striatella, and Anopheles gambiae mosquito.
Our work demonstrates that the combination of genomeresolved metagenomics, long-read sequencing, and pangenomic strategies provide an effective computational framework to investigate the diversity and distribution of mobile genetic elements in endosymbionts that are challenging to cultivate. Furthermore, it shows the importance of studying distinct individuals from wild mosquito populations, in parallel with controlled experiments in laboratory settings, to improve our understanding of the Wolbachia mobilome. The fragmented nature of Wolbachia MAGs in our study emphasizes the critical need for harnessing the power of emerging long-read sequencing technologies to characterize complex genomic variations of Wolbachia and its mobilome at finer scales. Wolbachia has been so far recalcitrant to`genetic modification, but the discovery of pWCP and phage WO may create new avenues for effective genome-editing strategies.

Methods
Sample collection. We collected mosquito specimens using a carbon dioxide mosquito trap located in Languedoc, Herault, France in May 2017 (Camping l'Europe de Vic La Gardiole, EID Méditerranée). Specimens were transported alive to the laboratory immediately upon recovery. We anesthetized adult females for 4 min at −20°C and proceeded to species-level identification. To remove potential contaminants from the insect surface, we gently vortexed specimens for 1 min in 1 ml cold (4°C) 96% ethanol. We then transferred them to a new 1.5 ml tube with 1 ml sterile cold (4°C) phosphate-buffered saline (PBS) 1× solution and gently vortexed them again for 10 s to avoid DNA precipitation with ethanol. Finally, we transferred specimens onto a sterile microscope slide with sterile PBS 1× on top of a cold plate and dissected two ovaries from four specimens using sterilized tweezers. We preserved ovary samples at −80°C until further processing.
C. pipiens complex controls. We obtained Wolbachia uninfected mosquitoes by treating the host with antibiotics. For antibiotics treatment we either used TC, which inhibits protein synthesis, or rifampicin, which interferes with nucleic acid synthesis 67 . The TC-treated lines 47 shown in Fig. 1 (C. quinquefasciatus SLAB-TC lines, ISEM, France) were raised at least 1 year without TC in standard laboratory conditions before beginning experiments. The rifampicin and oxytetracyclinetreated C. pipiens lines were kindly provided by Maria del Mar Fernandez de Marco (Animals Plant Health Agency, UK).
Metagenomic library preparation and sequencing. We extracted total genomic DNA from each ovary sample, hereafter referred to as O03, O07, O11, and O12, using the MoBio PowerFecal DNA Isolation Kit (QIAGEN Inc., Germantown, MD, USA). We used an E220 Covaris instrument (Covaris, Woburn, MA, USA) to sonicate 3.8-5.7 ng of genomic DNA. We end-repaired and 3'-adenylated resulting fragments and used the NEBNext Ultra II DNA Library Prep Kit for Illumina (New England Biolabs, Ipswich, MA, USA) to add NEXTflex PCR-free barcode adapters (Bioo Scientific, Austin, TX, USA). We purified ligation products by Ampure XP (Beckman Coulter, Brea, CA, USA) and PCR-amplified DNA fragments (>200 bp; 2 PCR reactions, 14 cycles) using Illumina adapter-specific primers and NEBNext Ultra II Q5 Master Mix (NEB). After library profile analysis using an Agilent 2100 Bioanalyzer (Agilent Technologies, Santa Clara, CA, USA) and quantitative PCR quantification using the KAPA Library Quantification Kit for Illumina Libraries (KapaBiosystems, Wilmington, MA, USA), we sequenced the library using a HiSeq4000 Illumina sequencer (Illumina, San Diego, CA, USA) at the Genoscope in Evry, France, generating 151 bp paired-end reads. To remove the least reliable data, we filtered the raw sequencing results using cluster intensity and chastity filter as described in ref. 68 .
Metagenomic assembly and binning. We used illumina-utils v1.4.4 69 to qualityfilter raw paired-end reads with the program 'iu-filter-quality-minoche' with default parameters, IDBA_UD v1.1.2 70 to assemble paired-end reads into contigs, and Bowtie2 v2.2.9 71 for all read recruitment analyses. We processed our contigs that are >1000 nts and read recruitment results with anvi'o v5 72 to recover MAGs from C. pipiens metagenomes. Briefly, we used the program 'anvi-gen-contigsdatabase' to generate anvi'o contigs databases for each four individual assemblies, during which anvi'o calculates and stores tetranucleotide frequency values and Prodigal v2.6.3 73 identifies ORFs in each contig. We ran default anvi'o HMM profiles on resulting contigs databases using the program 'anvi-run-hmms' and assigned functions to genes using the program 'anvi-run-ncbi-cogs', which searches gene amino acid sequences against the December 2014 release of the COG database 74 using blastp v2.3.0+ 75 . We profiled the short reads of our contigs recruited from each four individual metagenomes onto the four assemblies using the program 'anvi-profile', which generates anvi'o profile databases that store the coverage and detection statistics of each contig within each sample independently. We then merged resulting anvi'o profile databases for each sample using the program 'anvi-merge'. For an initial coarse binning, we used the CONCOCT 76 algorithm through the program 'anvi-cluster-with-concoct' and confined the number of clusters to five using the parameter '-num-clusters-requested 5'. We then used the program 'anvi-script-get-collection-info' to identify bins with a bacterial population genome, manually refined them to identify bacterial genomes using the program 'anvi-refine', and assigned taxonomy to resulting MAGs using NCBI's non-redundant protein sequence database with amino acid sequences of core genes. The program 'anvi-refine' allows the identification and refinement of population genome bins through an interactive interface by offering a number of tools including (1) guiding hierarchical clustering dendrograms of contigs (including one based on tetranucleotide frequencies and differential coverage statistics across samples), (2) real-time completion and redundancy estimates of binned contigs based on bacterial single-copy core genes 77 , (3) interfaces that display nucleotide-level coverage statistics per gene and contig, (4) functional annotations and synteny of genes, and (5) online sequence homology search options for gene and contig sequences. These tools collectively help minimize the inclusion of non-target contigs (such as eukaryotic contamination) in final bins, especially in complex metagenomes 78 . Except for the manual refinement step, we automated the processing of our raw sequencing data with a snakemake 79 workflow anvi'o v5 implements (http://merenlab.org/2018/07/09/anvio-snakemakeworkflows/), and the Code availablility section reports necessary configuration files to fully reproduce this analysis.
Sample preparation and long-read sequencing. To investigate the circularity of pWCP, we prepared for a long-read sequencing strategy. We extracted and pooled ARTICLE NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-019-08973-w DNA from the ovaries of 7 C. pipiens and 15 C. quinquefasciatus mosquitoes. As the total DNA quantity from this step was only 5.83 ng, we also extracted DNA from the abdomens of 15 C. quinquefasciatus that left us with a total of 76.63 ng. To acquire high-molecular weight DNA from these low-biomass samples, we used an in-house Phenol Chloroform protocol for total genomic DNA extraction. Briefly, we used a chloridric acid-washed pestle adapted to the size of a microcentrifuge tube to gently grind tissue in cetyltrimethylammonium bromide buffer (10 ml) mixed with betamercaptoethanol (20 µL) that we incubated for 15 min at 60°C. We then added 200 µL of phenol:chloroform:isoamyl alcohol and incubated each microcentrifuge tube at room temperature while rotation shaking for 5 min. To minimize inclusion of proteins, we centrifuged the tubes for 10 min at 8000 rpm (5939 × g), removed the aqueous phase (the top layer), and placed it into fresh tubes. We repeated those steps with chloroform:isoamyl alcohol and lastly with isopropanol. We incubated microcentrifuge tubes for 10 min at room temperature and mixed by inversion every 3 min. We then centrifuged those tubes for 10 min at 10,000 rpm (9279 × g), removed the top layer, and kept the sole DNA at the bottom of the microcentrifuge tube. We next washed DNA with 400 µL of freshly prepared 70% ethanol and centrifuged tubes for 10 min at 10,000 rpm (9279 × g). Finally, we kept microcentrifuge tubes in a speed vacuum for 10 min and eluted DNA with 20 uL of molecular grade water. During all laboratory steps, we performed smooth and slow pipetting with cut pipette tips to avoid further shearing the DNA molecules. To prepare mosquitoes' total DNA for sequencing, we used the Rapid Barcoding Kit (SQK-RBK004) by Oxford Nanopore Technologies. We followed the manufacturer protocols, except for one additional step, in which we supplemented samples that did not meet the manufacturer DNA input mass recommendations (400 ng) with linear double-stranded lambda DNA (New England Biolabs) for 'padding'. We then used 1× Agencourt AMPure XP beads (A63882, Beckman Coulter) for sample clean-up and concentration of pooled barcoded samples. A MinION with a pristine R9.4/FLO-MIN106 flow cell (Oxford Nanopore Technologies) sequenced the final prepared library with a starting voltage of −180 mV and a run time of approximately 20 h. While the flow cell had 1293 active pores prior to the run, only 200 were active at the beginning of the run, dropping down to 2 active pores by 20 h. We were not able to determine whether the low number of initial pores and their rapid decline was due to remaining chemical residues from our in-house Phenol Chloroform extraction protocol or a defective flow cell. We used MinKNOW (v1.15.4) for live base-calling and Albacore (v2.1.7) to de-multiplex our samples and convert raw FAST5 files to FASTQ files. For downstream analyses, we only used reads with a minimum quality score of 7.
Pangenomic analysis and the metapangenome. Our analysis of the Wolbachia pangenome followed the workflow outlined in ref. 80 and at http://merenlab.org/ 2016/11/08/pangenomics-v2/. Briefly, we first generated an anvi'o genomes storage database using the program 'anvi-gen-genomes-storage' for our four Wolbachia MAGs and the wPip reference genome wPip Pel 40 (NCBI Accession ID NC_010981.1). We then used the program 'anvi-pan-genome' using the flag '-usencbi-blast', and parameters '-minbit 0.5', and '-mcl-inflation 5'. This program (1) uses the blastp program all vs. all to create a graph of similarities between pairs of amino acid sequences, (2) removes bad hits using the 'minbit heuristic', (3) uses the Markov Cluster algorithm to determine gene clusters, (4) calculates the occurrence of gene clusters across genomes and the total number of genes they contain, (5) realizes a hierarchical clustering analysis for gene clusters (based on their distribution across genomes) and for genomes (based on shared gene clusters), and finally (6) generates an anvi'o pangenomic database that stores all results for downstream analyses. We used the program 'anvi-display-pan' to visualize gene clusters and their distribution across the five genomes in anvi'o interactive interface and summarized our manual selection of gene clusters using the program 'anvisummarize'. To identify 'phage-like' gene clusters, we used blastn to search for genes homologous to the five phage regions 27 in the wPip Pel genome (hereafter referred to as WOPip1-5). We finally used a custom R script (see the Code availability section) to recover coverage values of each gene cluster using our four metagenomes. We visualized the metapangenome using the program 'anvidisplay-pan'.
Computing the density of single-nucleotide variants per genome and metagenome. To infer the extent of heterogeneity in metagenomes for each Wolbachia MAG, we used the anvi'o program 'anvi-gen-variability-profile' to recover singlenucleotide variants from each merged profile without any filters (following the tutorial at http://merenlab.org/2015/07/20/analyzing-variability/). We then kept only single-nucleotide variants that fall into the context of complete gene calls to minimize erroneous variants that often emerge from mappig artefacts around beginnings and ends of contigs and calculated the percentage of nucleotide positions in each genome that was not clonal in the metagenome.
Computing the average nucleotide identity between genomes. To compute the level of similarity between the four Wolbachia MAGs and the wPip Pel reference genome, we used the program 'anvi-compute-ani', which called PyANI v0.2.7 81 in 'ANIb' mode to align 1 kbp long fragments of the input genome sequences with the NCBI's blastn program to summarize the average nucleotide identity and alignment coverage scores.
PCR amplification of DNA. We used a set of specific primers on the four C. pipiens individuals and three Wolbachia-free C. quinquefasciatus controls generated by TC treatment using standard techniques. The presence of arthropod DNA was verified by amplifying a ca. 708-bp fragment of the Cytochrome Oxidase I gene from arthropod mitochondria using LCO1490: (5'-GGT CAA CAA ATC ATA AAG ATA TTG G -3') and HCO2198: (5'-TAA ACT TCA GGG TGA CCA AAA AAT CA -3') primers 48 . We checked for the presence of Wolbachia DNA by PCR amplifying a ca. 438-bp fragment of the Wolbachia 16S rDNA gene using Wspec-F: CAT ACC TAT TCG AAG GGA TAG and Wspec-R:  AGC TTC GAG TGA AAC CAA TTC primers 49 . The PCR conditions for both sets of primers included a temperature regime of 2 min at 94°C, followed by 29 cycles of 30 s at 94°C, 45 s at 49°C, 1 min at 72°C, and a final elongation of 10 min at 72°C. To detect the circularity of the plasmid-like genome, we designed new primers (listed in Results) and used the following cycling parameters: 2 min at 95°C, followed by 34 cycles of (94°C for 30 s, 58°C for 45 s, and 72°C for 2 min) with a final elongation of 5 min at 72°C. Finally, to confirm the presence of the IS110 TE, we used primers described in the Results section under the following PCR conditions: 2 min at 95°C, followed by 34 cycles of (94°C for 30 s, 55°C for 45 s, and 72°C for 30 s) and a final elongation of 5 min at 72°C. Amplifications by PCR were performed using Platinium Taq DNA Polymerase High Fidelity (5 U/μL, Invitrogen).
Plasmid characterization. We used the Glimmer software within Geneious v8.1.9 (Biomatters, Ltd) to identify ORFs and used the NCBI databases for nonredundant protein sequences and conserved domains 82 , SMART 83 , and HHpred 84 (including SCOPe, Pfam, and COG) to manually annotate putative functions based on amino acid sequence homology searches. We used the IS finder platform (https://isfinder.biotoul.fr/) for IS homology search. The nucleotide analysis plugin in Geneious for EMBOSS identified and illustrated the EP in the plasmid. Code availability. The URL http://merenlab.org/data/2019_Reveillaud_and_ Bordenstein_et_al_Wolbachia gives access to a reproducible bioinformatics workflow document and ad hoc scripts used for all computational analyses.

Data availability
The raw sequencing data for shotgun metagenomes are available in the European Nucleotide Archive via accession code PRJEB26028. Sequencing data for the MinION run is available at https://doi.org/10.6084/m9.figshare.7306784. We also made available FASTA files for individual metagenomic assemblies (https://doi.org/10.6084/m9.