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 metagenomes42,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).

Table 1 Wolbachia MAG estimates Full size table

Wolbachia MAGs include highly covered non-phage contigs

The relatively low number of single-nucleotide variants (0.01–0.05%, Supplementary Table 2) suggested that each Wolbachia MAG represented a nearly monoclonal population of bacterial cells in the ovary metagenomes. In addition, the high average nucleotide identity across our MAGs and the 1482 kbp reference genome for Wolbachia, wPip Pel40 (99.1–99.98%, Supplementary Table 3) suggested a high degree of conservation between the endosymbionts of different individuals. Our metagenomic read recruitment analyses using Wolbachia MAGs revealed consistent coverage statistics averaging 168×–491× for contigs that were enriched with bacterial genes. However, a subset of contigs in each individual displayed approximately a five-fold increase in coverage (Supplementary Table 4). Because Wolbachia harbour prophage WO that can enter the lytic cycle and form phage particles27,44,45, we postulated that some of the contigs could be phage associated, which would explain coverage inconsistencies. We used the five prophage regions identified in the wPip Pel reference genome27,40 to identify contigs enriched with genes of phage origin. Contigs classified and validated through homology searches against phage WO matched to contigs that were highly covered in our MAGs, confirming that most shifts in coverage could be attributed to phages of Wolbachia. However, surprisingly, five contigs in our MAGs (Contig O12_A, Contig O11_A, Contig O07_A, Contig O07_A’, and Contig O03_A, Supplementary Data 1) showed no homology to prophage WO despite their remarkable coverage that ranged between 720× and 2176×. Interestingly, the 5’ and 3’ ends of these contigs showed homology to the non-coding flanking regions of wPip’s ISWpi12 transposable element (TE; WP0440, WP1209, and WP1347) of the IS110-family46. Given their (1) high coverage in metagenomes, (2) lack of homology to prophage WO, and (3) putative association with IS110 TE, we hypothesized that these contigs could represent extrachromosomal elements.

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 study32 confirmed the near identical presence of pWCP in all three, including the increase in coverage matching to the IS110 TE (Supplementary Figure 1).

Fig. 1 The artificially circularized genome of putative plasmid pWCP. a illustrates Contig O11_A and IS110 transposable element (TE) identified in our assembly. The red rectangles are regions of 100% nucleotide identity between the two contigs. Outward PCR primers were designed to amplify and confirm circularity of the sequence. b–e Gels for PCR tests to confirm a Wolbachia-associated circular genome. To verify the presence of arthropod DNA in our four Culex pipiens samples and the tetracycline-treated (TC) Culex quinquefasciatus samples, we PCR amplified a 708-bp sequence using LCO1490 and HCO2198 primers (b). A 438-bp fragment of the Wolbachia 16S rRNA gene (c), an approximately 1800 bp sequence amplified with the outward primers designed to support circularity of the genome, as illustrated in top panel (d) and a 431-bp of IS110 TE (e) were obtained in wild C. pipiens samples O03-O07-O11-O12 while no amplification was observed in Wolbachia-free samples. NC corresponds to negative control. f illustrates the complete genome. Each arrow represents an open reading frame (ORF). ORFs with no homology to a known function are shown in grey. ParA-like (green), RelBEtoxin–antitoxin operon (blue), and DnaB_C replicative DNA helicase (orange) that is disrupted by the ISWpi12 TE of the IS110 family (purple) are represented by arrows (with an e value < e−12 from an NCBI Conserved Domain or Pfam Search). Black squares represent the location of (1) the variable number tandem repeat (VNTR) and (2) the extragenic palindrome (EP) region Full size image

To validate the circular and extrachromosomal nature of pWCP independently of short-read recruitment- and assembly-based 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).

Fig. 2 Alignment of MinION long reads to pWCP. a The alignment of long reads that cover >50% of the pWCP genome (only 12 of 13 total long reads are shown; we omitted 1 from this display solely due to space considerations). Each rectangular figure shows high scoring pairs (HSPs) and their alignments between pWCP and long reads. The broken HSPs that are parallel to each other are due to low-quality regions in long read sequences, and they are shown in different shades. Concentric circles around pWCP demonstrate the alignment of each long read and their start and stop positions. b The alignment of long reads that cover <50% of the pWCP genome (only 5 of 6 are shown). Every long read shown in the figure has a hit to the IS110 TE Full size image

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 lines47. First, using LCO1490 and HCO2198 mitochondrial cytochrome c oxidase subunit I (COI) invertebrate primers48, 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 primers49. 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.150), 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 family46 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 genomes51 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 operon52,53. Most remaining pWCP genes were hypothetical and unique to either wPip and/or the B-Wolbachia phyletic supergroup54 (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 particles55. 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, b, Supplementary Figure 5), characterized as 16-nt 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 Wolbachia56,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 site56. 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 control58. Our analysis of the pWCP sequence also revealed a 209-bp extragenic palindrome (EP) region with two palindromes (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 palindromes34.

Fig. 3 pWCP contains a variable number tandem repeat (VNTR) region and extragenic palindrome (EP) sequence. a A VNTR region is located between parA and uncharacterized gene in the circular genome. While the number of repeats varies across individuals (b), the 36- and 33-bp spacers are conserved. Each black arrow represents a 16-nt repeat. The predicted DNA structure of the EP sequence is illustrated in c where color indicates probability of each base pairing. Red represents the strongest probability, whereas blue is the lowest Full size image

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.

Fig. 4 Wolbachia metapangenome in the context of wPip Pel genome synteny. The figure shows the presence–absence of 1166 gene clusters in the pangenome of four Wolbachia metagenome-assembled genomes (MAGs) and the reference genome wPip Pel. The gene clusters (i.e., groups of homologous genes based on amino acid sequence identity) are organized based on wPip Pel synteny. Each MAG is represented by two layers, where the first layer indicates the presence or absence of a gene cluster in a given MAG, and the other shows the average coverage of each Wolbachia MAG O07 gene cluster in the corresponding C. pipiens metagenome. The second to last layer shows whether genes in a given gene cluster have a match in NCBI’s COGs. The outermost layer associates gene clusters with previously identified prophage regions in the wPip Pel genome. The number of gene clusters assigned to WOPip prophage regions are indicated in parenthesis Full size image

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 wasps27. 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 pathogens59,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 associated50,61,62,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 chromosome40. 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 closed40, 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 functions30. 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 WO27.

Metagenomic read recruitment revealed between a 1.5- and 5-fold 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 mosquitoes44,45, observed lytic activity of Wolbachia phages in Nasonia testes, and the sequencing of WO genomes from purified phage particles27. 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 T464,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 genomes61,66 as well as the presence of TEs within and/or flanking prophage variants27,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.