Keywords:Wolbachia, evolution, endosymbiont, genomics

Cell Lines and Maintenance ofWolbachia in Cell Culture

This study used two mosquito cell lines that were infected by theWolbachia strainwMelPop and maintained with continuous passaging for several years (McMeniman et al. 2008). Briefly,wMelPop bacteria were purified fromDrosophila melanogaster w¹¹¹⁸ embryos (Min and Benzer 1997;Dobson et al. 2002) and established in aWolbachia-freeA.albopictus mosquito cell line, Aa23-T (O’Neill et al. 1997).wMelPop was maintained in Aa23 cells for a period of approximately 27 months before purification and subsequent introduction into a second mosquito cell line, RML12, also derived fromA.albopictus (McMeniman et al. 2008). ThewMelPop-infected RML12 cell lines were maintained for a period of 17 months prior to purification and DNA extraction for genome sequencing. Throughout the 44 months thatwMelPop was maintained in cell culture, insect cells were passaged every 3–4 days and approximately 20% of the cells were used to establish the next generation. Throughout this time, aliquots ofWolbachia-infected cells were collected and stored in liquid nitrogen.

DNA Purification ofwMelPop-CLA from Cell Culture and Mosquitoes

ForwMelCS,wMelPop, andwMelPop-PGYP, purification procedures were followed as described inIturbe-Ormaetxe et al. (2011). Note thatwMelPop-PGYP is referred to in that paper as “wMelPop-CLA fromA. aegypti PGYP1.”wMelPop-CLA was purified from cell lines as follows. To obtain enough material for the purification ofWolbachia DNA from cell lines, 20–30 175 cm² flasks containing confluent monolayers of cells were harvested after gently shaking the bottles. Cells were centrifuged in 50-ml conical flasks at 3,200 × g for 10 min at 4 °C. Culture media was discarded and the pellets were washed twice by resuspending them in SPG buffer (0.25 M sucrose, 0.2% BSA, and 10 mM MOPS, pH 7.2) and then centrifuged at 3,200 × g for 10 min at 4 °C. The supernatant containingWolbachia was sequentially filtered through 5-, 2.7-, and 1.2-µm syringe filters, pelleted by centrifugation at 18,000 × g for 20 min, and then resuspended in cold SPG buffer. Six hundred microliters of resuspendedWolbachia were carefully layered on top of a discontinuous Percoll gradient composed of 2.7 ml of 10% (v/v), 3 ml of 15%, 2 ml of 35%, and 4 ml of 50% Percoll/SPG (Duplouy et al. 2013). The equivalent of six to eight bottles of cells was used for each gradient tube. Tubes were centrifuged for 1 h at 8,700 rpm in a Beckman Optima-L-80 XP ultracentrifuge using a swinging bucket rotor SW41. Following centrifugation, four bands containing white cellular material were clearly visible in the interphase between the Percoll layers. The material in each band was recovered by removing the liquid above and pipetting the band out sequentially without disrupting the gradient.

For the extraction of DNA from the different bands, 750 µl of cell material was transferred to Eppendorf tubes and treated with DNaseI for 30 min at 37 °C to remove host DNA contamination. After treating the bands with 1 µl RNase-DNase free (Fermentas),Wolbachia cells were disrupted by incubation with proteinase K at 56 °C for 30 min. DNA was extracted using phenol/chloroform, precipitated, washed in 70% ethanol, and resuspended in TE or milli Q water (Millipore).

Total DNA was quantified using a nanodrop spectrophotometer. An aliquot containing approximately 500 ng of the obtained DNA was run on an agarose gel to test DNA quality and purity. DNA that was isolated from each of the different bands in the gradient was characterized for the presence of host, mitochondrial and bacterial contamination by polymerase chain reaction (PCR) using specific primers listed insupplementary table S1 (Supplementary Material online). Band 4 contained the highest concentration ofWolbachia, very low levels of host nuclear DNA, and no detectable mitochondrial contamination, and this band was therefore chosen for genomic DNA extraction and sequencing.

Genome Sequencing, Assembly, and Annotation

Initial genome sequencing was performed by the Australian Genome Research Facility. ThewMelCS genome was sequenced using Illumina, whilewMelPop,wMelPop-CLA, andwMelPop-PYGP were sequenced using 454 (Roche) pyrosequencing (seetable 1 for details). Subsequent Illumina resequencing ofwMelPop andwMelPop-PGYP was performed by the Ramaciotti Centre, University of New South Wales.

ThewMelPop,wMelPop-CLA, andwMelPop-PGYP genomes were each assembled using Newbler v2.6. The initial assemblies contained a substantial number of homopolymer errors, which we corrected forwMelPop andwMelPop-PGYP using Illumina sequencing data generated from the same DNA material that was used for the original 454 sequencing. For each genome assembly, we mapped to the assembly the 454 and Illumina reads from that strain. We then used Nesoni (http://www.vicbioinformatics.com/software.nesoni.shtml, last accessed November 19, 2013) to call variants for each mapping. We considered a variant to be evidence of a homopolymeric sequencing error in the assembly if the following two conditions were met: 1) the Illumina reads, which should not be subject to systematic homopolymer errors, were consistent with the variant and inconsistent with the assembly, and 2) there was disagreement about homopolymer length in the 454 reads mapped to the variant site. In these cases, we corrected the variant to match the Illumina data. Once assembly correction was complete, we compared the genomic arrangement of each strain with that of the completewMel genome using Mauve (Darling et al. 2004).

The correctedwMelPop genome assembly was automatically annotated using the NCBI Prokaryotic Genome Annotation Pipeline (PGAAP;http://www.ncbi.nlm.nih.gov/genomes/static/Pipeline.html, last accessed November 19, 2013). We then compared the gene complements ofwMel andwMelPop by performing reciprocal BlastN analyses of the genes annotated in the two strains.

Phylogenetic Analysis

To construct a whole-genome phylogeny ofwMel,wMelCS, andwMelPop, we used the correctedwMelPop genome assembly, the publishedwMel genome (Wu et al. 2004), eight of thewMel consensus genomes and the twowMelCS consensus genomes called byRichardson et al. (2012), and a consensus genome generated from ourwMelCS/Canton-S sequencing data using the same method asRichardson et al. (2012). Briefly, we mapped thewMelCS reads to the referencewMel genome using BWA (Li and Durbin 2009), then generated a pileup file with minimum and maximum read depths set to 10 and 100, respectively, and converted the resulting fastq file to fasta format. We then used Mauve (Darling et al. 2004) to align thewMel andwMelCS genomes to thewMelPop assembly, and exported the core genome regions with minimum LCB length of 100. This produced an alignment of 1,136,727 nt. We then inferred a maximum likelihood phylogenetic tree using RAxML (Stamatakis 2006), with a general time reversible model of nucleotide substitution with a gamma model of rate heterogeneity with four rate categories.

Identification of Sequence Variants between Genomes

We identified sequence differences between each of the genomes analyzed here (wMel, the threewMelCS genomes,wMelPop,wMelPop-CLA, andwMelPop-PGYP) using three main techniques:

1. For the threepopcorn genome assemblies, we used Mauve (Darling et al. 2004) to create pairwise alignments of each draft assembly to the other assemblies and to thewMel genome, and then exported the core alignments using the stripSubsetLCBs script (provided by Mauve developers athttp://gel.ahabs.wisc.edu/mauve/snapshots/, last accessed November 19, 2013). We then used custom Perl scripts to identify mismatches in the alignments.

2. To identify variants between the threepopcorn genomes, and between the threewMelCS genomes and thepopcorn genomes, we mapped the reads of each of the six data sets to each of thepopcorn assemblies, and called variants as described later.

3. We also performed an indirect comparison between these six data sets by mapping the reads of each data set to thewMel genome, calling variants, and then comparing variant calls across strains.

For the second and third approaches, we called variants using several methods depending on the kind of sequence data available. For 454 data sets, we mapped reads using Newbler and then examined high confidence variant calls. Illumina data sets were mapped and variants called using two complementary approaches, as follows.

First, reads were aligned with BWA (Li and Durbin 2009) usingaln andsampe with default parameters, and duplicates were removed withrmdup. SAMtools (Li et al. 2009)mpileup (with parameters -C50 -BEA) andbcftools (with parameter -D220) were run to call variants and produce VCF output. Variant sites were then filtered for minimum quality of at least ten and read depth between 10 and 220. Variant calling was repeated on the BWA mapping using Freebayes (Garrison and Marth 2012), and the same filtering steps were applied.

Second, variants were called with the Nesoni high-throughput sequencing data analysis toolkit (http://www.vicbioinformatics.com/software.nesoni.shtml, last accessed November 19, 2013), using SHRiMP (David et al. 2011) as the aligner and Freebayes as the caller.Nesoni samshrimp was run with default parameters, andnesoni filter was run (with parameter −−monogamous no to retain reads mapping to more than one location). Variants were called by first runningnesoni freebayes with parameters −−depth-limit 220 and −−ploidy 4 and then reducing the ploidy to 1 by runningnesoni vcf-filter on the VCF file produced by Freebayes.

Variants identified in one strain were also checked in all other strains, both bioinformatically (using the pipelines above and via manual inspection of read alignments) and using Sanger sequencing. By doing this, we hoped to eliminate biases in variant detection caused by differences in sequence quality or depth in different strains.

Copy Number Variation across the Genome

We searched for large-scale variations in copy number across the genomes by mapping the sequencing reads of each data set to thewMel genome using BWA (for Illumina data) or Newbler (for 454 data). For the coverage plots infigure 3, we calculated mean per-site total coverage (i.e., including assignment of multiply-mapping reads to a random instance of the repeat in the genome) for nonoverlapping 50-nt windows along the genome. Once we had identified the region of coverage variation in these genomes, we then used additional analyses to identify the boundaries of the triplication and deletion more precisely, as described insupplementary information (Supplementary Material online).

To confirm the triplication of the 19-kb region inwMelPop and its deletion inwMelPop-PGYP, we determined the relative PCR amplification of two unique sequences in that region (sequences spanning gene boundaries WD0512-WD0513 and WD0513-WD0514) compared with the single copywsp gene as a reference. As a control, we also calculated the relative PCR amplification of the single copy gene WD1213 compared towsp. DNA was extracted fromwMelCS- andwMelPop-infected flies andwMelPop-PGYP-infected mosquitoes using a DNeasy Blood and Tissue kit (Qiagen), and 10 ng of DNA was used for qPCR amplification using the primers described insupplementary table S1 (Supplementary Material online). The qPCR reaction contained 10 ng DNA, 5 µl of 2X LightCycler 480 Probes Master containing SYBR green (Roche), and 1 µM of each primer in a total volume of 10 µl. Reactions were performed in triplicate in a LightCycler 480 Instrument (Roche) with the following conditions: 95 °C for 5 min, and 45 cycles of 95 °C for 10 s, 60 °C for 15 s, and 72 °C for 1 s. Relative amplification was calculated using the formula E^(Cpwsp)/E^(Cp test gene), where E is the qPCR amplification efficiency, Cp the crossing point, and the test gene was WD0512–WD0513, WD0513–WD0514, or WD1213.

Searching for Other Structural Differences or New Genes

We used a number of approaches to search for structural variants and insertion of novel genetic material in these genomes. First, we inspected the high confidence structural rearrangements predicted by Newbler after mapping each of thewMelPop sequencing reads to thewMel genome. Four variants were identified, all of which were already known: the large inversion, the large triplication, duplication of ankyrin repeats in WD0766, and insertion of an IS5 element in WD1310. Second, we searched specifically for evidence of any additional movement of IS5 elements inwMelCS or thepopcorn genomes. We used the IS5 terminal inverted repeat sequence as a blastN query against the sequencing reads from each data set, and identified reads that matched both this repeat sequence and flanking unique sequence. These reads were then mapped to thewMel genome to identify insertion sites of all IS5 elements from each strain.

Finally, we mapped the reads forwMelCS,wMelPop andwMelPop-PGYP to thewMel genome, then used Perl scripts to identify three sets of cases: 1) read pairs that mapped significantly further apart than expected, possibly due to either genomic rearrangement or a deletion in the query genome, 2) read pairs that mapped in an unexpected orientation, possibly indicating genomic rearrangement, and 3) read pairs in which one read mapped towMel and the other did not, which could be due to insertion of novel genetic material into the query genome. For each set of cases, we mapped the reads in these categories to thewMel genome and identified regions where there was a concentration of mapped reads. Alignments at all potential sites of interest were then inspected manually. No additional structural variants or sites of insertion of novel genetic material were confirmed.

Timing of Mutations

To determine the time of appearance of the five mutation events detected during cell passage, cell stocks frozen in liquid nitrogen at different times during cell culture were revived and DNA extracted for PCR analysis. The presence of the 57-kb deletion was detected by PCR amplification of four single copy genes in that fragment (WD0511, WD512, WD513, and WD0514). The insertion of an IS5 element between WD0765 and WD0766 was determined by PCR using primers flanking the insertion site, whereas the presence of the 10 bp in WD0413 was confirmed using a PCR primer whose 3′ end sits in the deletion and only amplifies from the wild type sequence (primers listed insupplementary table S1,Supplementary Material online). Following the screening of all the frozen stocks, the PCR bands obtained were confirmed by sequencing (fig. 6).

Draft Popcorn Genomes

The draftwMelPop genome consists of 13 scaffolds, ranging in size from 2,302 to 547,521 nt, and an additional 12 unscaffolded small contigs between 599 and 1,946 nt (table 1). The draft genome ofwMelPop-PGYP is of similar completeness, consisting of ten scaffolds (2,300 to 541,636 nt in length) and 16 unscaffolded contigs (503–1,462 nt). ThewMelPop-CLA genome was sequenced using shotgun rather than paired-end reads, and the assembly consists of 220 unscaffolded contigs, between 504 and 34,297 nt in length. Like otherWolbachia genomes, thepopcorn genomes are strongly AT-biased, with an average of 36% GC content.

We annotated only thewMelPop draft genome. The annotation contains 1,111 protein-coding genes. All genes annotated inwMel have orthologs (or paralogs collapsed into a single contig, in the case of repeat genes such as IS5 elements) in thewMelPop assembly. A comparison of thewMelPop assembly with thewMel genome confirmed three previously described genomic differences between them. First, although the genomes of bothwMel andwMelPop each have 13 copies of the IS5 mobile element, only 12 of these copies are shared by both strains. An IS5 element present inwMel between genes WD0516 and WD0517 is not present inwMelPop; conversely, an IS5 element has been inserted inwMelPop into the ortholog of thewMel gene WD1310 (Riegler et al. 2005;Cook et al. 2008). Second, ankyrin repeat domains have been duplicated in the orthologs of WD0550 and WD0766 inwMelPop, and there are also differences in repeat number in the variable tandem repeat region VNTR-141 between the two strains (Riegler et al. 2005,2012). Third, a ∼143 kb region of the genome has been inverted betweenwMel andwMelPop (Sun et al. 2003;Riegler et al. 2005). The inversion encompasses genes orthologous to WD0399 to WD0536, and is flanked by identical inverted repeat sequences, atwMel coordinates 376,721–379,499 and 522,494–525,272. We also identified over 150 novel differences betweenwMel andwMelPop. All variants discussed later were checked using Sanger sequencing inwMel,wMelCS,wMelPop,wMelPop-CLA, andwMelPop-PGYP strains.

We identified and confirmed 156 single nucleotide changes or small indels between the annotatedwMel genome sequence and the draftwMelPop genome (supplementary data table S1,Supplementary Material online). Of these, 112 (104 SNPs and 8 indels) occur within putative coding regions, and 78 of the SNPs result in a coding change. Five frameshift changes were identified. The first, in the ortholog of WD1155, which encodes a hypothetical protein, should result in the production of a slightly longer protein (118 aa inwMelPop vs. 101 aa inwMel). Each of the remaining four frameshifts results in the loss of a stop codon, leading to the production of a single gene inwMelPop from what are annotated as two contiguous genes inwMel. In each case, thewMelPop CDS has full-length blast matches to genes in otherWolbachia strains, suggesting that the single-CDS state is ancestral, and these genes have been relatively recently pseudogenized inwMel. These four genes are the orthologs of 1) WD0026–WD0027, encoding a hypothetical protein, 2) WD1043–WD1044, a second hypothetical protein, 3) WD1215–WD1216, encoding a sensor histidine kinase/response regulator, and 4) WD1231–WD1232, encoding the protoheme biosynthesis proteinHemY.

Relationship betweenwMel,wMelCS, andwMelPop

Earlier analyses based on small numbers of sequence differences suggested thatwMel andwMelPop were sister strains andwMelCS a slightly more distant relative (Paraskevopoulos et al. 2006), while analysis of the large inversion, IS5 insertion sites and VNTRs indicated thatwMelCS andwMelPop were more closely related to each other than towMel (Riegler et al. 2005,2012). To gain a fuller understanding of the relationship between these strains, we wished to compare their complete genomes. We did not have an assembled genome sequence ofwMelCS, but we used Illumina sequence data from threewMelCS lines. We purified the first of these from a lab line ofD. melanogaster Canton-S flies (Iturbe-Ormaetxe et al. 2011). The second and third were identified byRichardson et al. (2012) from the sequencing reads of twoD. melanogaster lines collected from a single population in Raleigh, North Carolina, in 2003 and sequenced as part of theDrosophila melanogaster Genetic Reference Panel (Mackay et al. 2012).Richardson et al. (2012) generated a consensus genome sequence for each of these twowMelCS lines (and for many additionalwMel lines) by mapping the sequence reads from each line to thewMel genome; for consistency we called a consensus genome sequence of ourwMelCS line using the same technique. We then aligned the publishedwMel genome, twowMel consensus genomes from each of the fourwMel clades identified byRichardson et al. (2012), the threewMelCS consensus genomes, and the draftwMelPop genome assembly, and constructed a maximum likelihood phylogenetic tree of the sequences (fig. 2).

wMelPop clusters clearly withwMelCS to the exclusion of allwMel sequences.wMelPop andwMelCS from Canton-S flies (shown aswMelCS/Canton-S) are extremely closely related; the two Raleigh sequences (shown aswMelCS/DGRP335 andwMelCS/DGRP338) cluster together on a separate branch. This branching pattern indicates thatwMelPop is recently derived from withinwMelCS. These inferred relationships are not an artifact of generating thewMelCS consensus genome sequences by mapping to thewMel genome: The same patterns of relatedness are observed whenwMelCS reads are mapped to thewMelPop assembly for the purposes of variant calling.

Genomic Basis of Life-Shortening: Differences betweenwMelCS andwMelPop

wMelCS is the nonpathogenic strain most closely related towMelPop. This means that the genomic changes that causedwMelPop to become pathogenic must have occurred betweenwMelCS andwMelPop. To attempt to identify these changes, we mapped the three sets ofwMelCS reads to thewMelPop draft assembly and searched for differences between them. We found a small number of variants between thewMelCS/DGRP335 andwMelCS/DGRP338 sequences andwMelPop (table 2), but at these sites,wMelCS/Canton-S always matched the sequence ofwMelPop, indicating that these differences could not be associated with pathogenesis.

ThewMelPop draft genome assembly contains collapsed repeats and could possibly contain other errors that might impede the detection of sequence variants. To ensure that we were not missing true variants, we therefore also mapped the sequencing reads from each of the threewMelCS data sets and thewMelPop,wMelPop-CLA andwMelPop-PGYP data sets against thewMel genome, and compared the variants that were called. We again identified differences between the DGRPwMelCS sequences andwMelPop (table 2), but no SNPs or indels that could differentiatewMelCS/Canton-S fromwMelPop. We also found no evidence thatwMelPop andwMelCS differ in the insertion sites of IS5 transposable elements or other mobile elements.

We did, however, detect a region of copy number variation betweenwMelCS/Canton-S andwMelPop. A region of thewMelPop genome corresponding to the genes WD0507 to WD0514 inwMel has sequence coverage approximately three times higher than that of the rest of the genome (fig. 3;supplementary fig. S1,Supplementary Material online). There is no obvious variation in coverage of any of thewMelCS genomes in this approximately 19-kb region. To confirm the triplication of these genes inwMelPop, we performed qPCR experiments comparing normalized amplification of genomic DNA fromwMelCS in Canton-S flies,wMelPop inw¹¹¹⁸ flies andwMelPop-PGYP inA. aegypti. Using primers that spanned the gene boundaries WD0512–WD0513 and WD0513–WD0514, and normalizing against the single-copy genewsp, we found that these genes amplify approximately three times as highly inwMelPop as inwMelCS (fig. 4). These results are consistent with earlier data serendipitously examining expression of these genes using Southern blots (Iturbe-Ormaetxe et al. 2005): in those data, the hybridization signal for these genes, and in particular for WD0514, is stronger inwMelPop than inwMel orwMelCS, whereas the other genes tested, includingwsp, produced similar signal intensities for the different strains. It seems unlikely, however, that the triplication of one or all of these genes could be directly responsible for pathogenesis, as the region that is triplicated inwMelPop is completely deleted inwMelPop-CLA andwMelPop-PGYP (confirmed by mapping of sequencing reads to thewMel genome for both substrains,fig. 3 andsupplementary fig. S5,Supplementary Material online, and with qPCR forwMelPop-PGYP,fig. 4), and yet these two substrains remain pathogenic.

We also tested whether the genetic basis of pathogenesis might be associated with a plasmid present inwMelPop but not the nonpathogenic strains. If such a plasmid were present, we would expect thewMelPop assembly to contain one or more contigs that do not correspond to matching sequence inwMel. No such contigs were identified in the assembly. We also attempted to separately assemble all reads from thewMelPop sequencing run that did not map to thewMel genome, but this produced no large contigs that might form part of a novel plasmid. This is in agreement with earlier laboratory work bySun et al. (2001) who found no evidence of extrachromosomal DNA on pulse field gel electrophoresis (PFGE) ofWolbachia strains includingwMelPop, indicating that plasmids were not present.

Given the very limited genomic differences we have detected betweenwMelCS andwMelPop, we also wished to confirm that the phenotypic differences betweenwMelCS andwMelPop are due solely toWolbachia, and not to host effects.wMelPop retains its pathogenic phenotype even after being purified and microinjected intoA. aegypti (McMeniman et al. 2009), demonstrating that the pathogenic effects associated with this strain are clearly caused byWolbachia and are not due to host nuclear or mitochondrial factors. However, it is possible, if unlikely, thatwMelCS might also have pathogenic capabilities that are suppressed in some way by the Canton-S host background. To test this, we purifiedwMelPop and transinfected it into Canton-S flies that had been antibiotic-treated to remove their nativewMelCS infection. If Canton-S flies are capable of suppressing pathogenesis, then these flies should live as long as Canton-S flies carryingwMelCS. But they do not: female Canton-S flies infected withwMelPop (at generation G40 after transinfection) die significantly earlier than Canton-S flies infected withwMelCS (median survival 26 days vs. 55.5 days; Cox regression,χ² = 4.72, df = 1,P = 0.030; seesupplementary information [Supplementary Material online] for further details). This confirms that the difference in pathogenicity betweenwMelCS andwMelPop is due not to host effects but toWolbachia strain, despite the few genomic differences we have observed between these strains.

Rapid Evolution ofwMelPop after Transinfection into Cell Lines

As part of a strategy to attempt to preadaptwMelPop derived fromD. melanogaster before establishing an infection in the mosquitoA. aegypti,wMelPop was passaged in mosquito cell lines for approximately 44 months before being reintroduced to flies and phenotypically recharacterized (McMeniman et al. 2008) (fig. 1). We compared the draft genomes ofwMelPop and the cell-line passaged strainwMelPop-CLA to identify the genetic changes that occurred during evolution in a novel cellular environment, which may be associated with the subsequent attenuation of pathogenesis inD. melanogaster. Five genetic differences were detected: an IS5 insertion, a multigene deletion, two point mutations, and a 10-bp deletion.

Timing of Changes in Cell Lines

During serial passaging, mosquito cells infected withwMelPop were passaged every 3–4 days; cells were periodically snap frozen and stored in liquid nitrogen as part of the usual maintenance routine for tissue culture. This collection of frozen cells provided an opportunity to estimate when genetic changes occurred during the evolution ofwMelPop in mosquito cells. Unique sets of PCR primers (supplementary table S1,Supplementary Material online) were used to screen these banked cells for the mutations discovered in thewMelPop-CLA genome.

The first mutation to occur was the IS5 element insertion between the orthologs of WD0765 and WD0766, which was detected only 13 months afterwMelPop was established in the Aa23 cell line (fig. 6). For a period of 9 months, twowMelPop variants could be detected in the Aa23 cell lines: the wild type (IS5 absent from locus) or thewMelPop-CLA form (IS5 present at locus). By 21 months after infection of Aa23, the IS5 insertion at this locus had become fixed within theWolbachia population.

The 57-kb deletion was detectable 25 months after transinfection and was fixed within theWolbachia population 8 months later; during this period,Wolbachia was purified from the Aa23 cell line and introduced into a secondA. albopictus cell line, RML12 (fig. 1). As both the wild type and mutant forms were detected in the early RML12 cell cultures, it is unlikely that the deletion became fixed as a result of a population bottleneck imposed due toWolbachia transinfection into a new insect cell line. Finally, the 10-bp deletion within WD0413 was detectable 37 months after initial cell line infection (10 months after transinfection form Aa23 to RML12), but only became fixed within the population 15 months later.

Evolution ofwMelPop-CLA after Transinfection into Mosquitoes

To test whether this rapid rate of genomic change continued afterwMelPop-CLA was transinfected into mosquitoes, we also purified and sequenced the genome of thewMelPop-PGYP strain from mosquitoes (McMeniman et al. 2009) 4 years after the infection was introduced into this host. We compared thewMelPop-CLA andwMelPop-PGYP genomes by 1) aligning the assemblies to one another and identifying mismatches, 2) mapping the reads of each strain against the assembly of the other strain, and 3) mapping reads of both strains to thewMel genome, then calling and comparing variants. We observed no genetic differences betweenwMelPop-CLA andwMelPop-PGYP. We also found no evidence of novel polymorphisms segregating in the population ofwMelPop-PGYP sequenced.

Origin ofwMelPop

The Canton-S line ofD. melanogaster was collected in Canton, Ohio, prior to 1938 (Bridges and Brehme 1944), and has since been maintained as a common laboratory stock. ThewMelCSWolbachia strain carried by Canton-S has at least eight genomic differences (table 2) from thewMelCS found in the DGRP335 and DGRP338 lines that were collected in 2003 in Raleigh, North Carolina (Mackay et al. 2012). In contrast, we identified only one difference (a change in copy number of a single genomic region) between Canton-SwMelCS and the pathogenicwMelPop. The original location and date of collection of the line carryingwMelPop is unclear, but it had been established in the laboratory prior to 1948 (Hannah 1949;Valencia and Muller 1949). The inferred collection dates ofwMelPop andwMelCS/Canton-S and the extremely close genomic similarity between these two strains suggest that the pathogenicwMelPop strain arose from within thewMelCS clade at some time in the mid-20th century.

It is possible that the mutation/s that led towMelPop becoming pathogenic occurred in the wild before the line carrying this strain was collected. However, as the fitness costs ofwMelPop are high, it seems more likely that the evolution of pathogenesis occurred in the laboratory after collection of a line carrying a benignwMelCS strain. Although theD. melanogaster line carrying the proto-wMelPop was crossed with irradiated males of other lines (Hannah 1949;Valencia and Muller 1949), there is no evidence that females of this line were directly exposed to mutagenizing agents. In the absence of paternal inheritance ofWolbachia, the mutation/s that led to the development of pathogenesis are likely to have arisen as the result of normal errors in genomic replication, and to have been maintained due to relaxed selection for longevity in a fly stock center environment.

Genomic Basis of Pathogenesis ofwMelPop

Despite their dramatic differences in phenotype, we identified only a single genomic difference betweenwMelCS andwMelPop: the triplication in copy number of a 19-kb genomic region. This region, which has previously been shown to be highly labile (Iturbe-Ormaetxe et al. 2005;Riegler et al. 2005;Woolfit et al. 2009), contains eight genes, most of which are either transposon-related or annotated as hypothetical proteins. However, it seems unlikely that this increase in copy number is itself associated with pathogenesis, as this same genomic region has been deleted from the pathogenic substrainswMelPop-CLA andwMelPop-PGYP.

Do the genome sequences ofwMelCS andwMelPop strains differ in other ways that we have not discovered? We know that our data have sufficient power to identify many sequence differences: we found over 150 single nucleotide changes and indels that varied betweenwMel andwMelCS/wMelPop. More than 90% of these were independently called as high-confidence variants using sequence data from each ofwMelCS andwMelPop, and all variants that were called in only one data set could later be identified in the sequence reads of the other data set. We were also able to draw on data from three independent sequence data sets for bothwMelCS andwMelPop strains, making it unlikely that stochastic variation in sequence coverage for one genome might be obscuring a genomic difference between them. Nonetheless, there are a number of types of sequence variants that could have remained undetected by our analyses.

First, the genomes ofwMel,wMelCS, andwMelPop are rich in sequence repeats, from short tandem repeat sequences to multiple copies of transposase-related genes each over 2,000-nt long (Wu et al. 2004;Cerveau et al. 2011;Leclercq et al. 2011;Riegler et al. 2012). These repeats present difficulties for both read alignment and de novo assembly (Treangen and Salzberg 2012), and it is possible that sequence variants present in a subset of repeat copies inwMelCS and/orwMelPop have not been identified. The great majority ofwMel annotated repeat sequences longer than 200 nt are associated with mobile elements (Wu et al. 2004), and it seems functionally unlikely that minor sequence variants in one or more copies of these genes might cause pathogenesis. Shorter repeat sequences, however, have been linked to multiple modes of pathogenesis (Delihas 2011;Treangen and Salzberg 2012), and undetected variation in similar sequences inwMelPop could be contributing to its pathogenic phenotype.

Second, the fact that numerousWolbachia repeat sequences are longer than the insert size of our paired-end reads means that it is not always possible to assemble or align reads across repeat regions. If long repeats form the breakpoints of structural rearrangements, these events might not be detected. Two observations argue against the possibility that we have overlooked any large genomic rearrangements, however: 1) A large inversion present in thewMelPop genome was successfully identified in our mapping analyses, despite being flanked by repeats of moderate length, and 2) an earlier physical and genetic map ofwMelPop based on restriction endonuclease digestion identified that same inversion as the only large-scale disruption of collinearity betweenwMelPop andwMel (Sun et al. 2003).

Third, although we identified a number of indels in our data sets, indel detection is more challenging and currently less accurate than SNP calling for next-generation sequence data (Albers et al. 2011), and additional indels may have been missed. An indel resulting in a frameshift within a coding gene, or disrupting an intergenic regulatory region, could lead to changes in protein function or expression. Transcriptomic analyses would provide a more direct and powerful way of detecting such changes. Finally, next-generation sequencing read coverage of genomes is nonrandom but shows biases associated with GC content, proximity to the origin of replication and other factors (Minoche et al. 2011;Meglecz et al. 2012), and so someWolbachia genomic regions may be systematically underrepresented in our data. However, when we map thewMelCS orwMelPop reads against thewMel genome, more than 99% of bases in the reference have better than 10X coverage, which should provide sufficient data to detect variants if present.

It is possible thatwMelCS andwMelPop differ more at the epigenetic than the genomic level. Previous work has shown that bacteria possess diverse adenine methylation systems with a wide range of specificities and activities (Low et al. 2001;Murray et al. 2012), and that changes in this methylation can have large-scale effects on bacterial gene regulation (Fang et al. 2012) and modulate bacterial phenotypes including virulence (Low et al. 2001;Heusipp et al. 2007). The genomes of numerousWolbachia strains are known to encode two phage-derived adenine methylases (Saridaki et al. 2011), and homologs of these genes are present inwMel,wMelCS, andwMelPop, suggesting that these strains have the genetic machinery required to differentially methylate their genomes. Genome-wide analyses of methylation are becoming increasingly tractable (Murray et al. 2012), and this is a promising avenue for future research on the differences betweenwMelCS andwMelPop.

Regardless of potential epigenetic differences betweenwMelCS andwMelPop, however, our results demonstrate that endosymbiotic bacteria that differ very little at the genomic sequence level can cause extremely different phenotypes in their hosts.

Rate of Genomic Change after Transfer to a New Host

wMelPop and thewMelCS strain from Canton-S flies have been evolving separately for at least 70 years, and yet show very little detectable genetic divergence. This is in strong contrast to the rapid evolution we observed inwMelPop after it was transferred into mosquito-derived cell lines. The first mutation, the movement of an IS5 element, occurred within 13 months of transinfection into cells; in contrast, we have identified only two changes in IS5 element location betweenwMel andwMelCS, which diverged several thousand years ago (Richardson et al. 2012). Subsequent genomic deletions and single nucleotide changes occurred within 4 years of the initial transfer to cell lines. This rapid evolution may be due to adaptation to a new host (fly to mosquito), adaptation to the cell line environment, and/or drift due to relaxed selection in cell lines.

There have been relatively few previous studies comparing the complete genome sequences of bacteria before and after transfer to a new host, and most have examined strains or species that have diverged for much longer time periods than the 4 years that separatewMelPop andwMelPop-CLA. Nonetheless, a number of common patterns associated with bacterial host jumps have emerged. Two independent transfers ofStaphylococcus aureus from humans to novel hosts have occurred recently: to ruminants within the last 115–1,204 years (Guinane et al. 2010) and to poultry 30–63 years ago (Lowder et al. 2009). In both cases, host adaptation has occurred via a combination of gene loss, acquisition of horizontally transferred genes, and a small-to-moderate amount of diversification in gene sequences. An older transfer, ofHelicobacter from humans to large felines estimated to have occurred 50,000 to 400,000 years ago, shows the same pattern: host adaptation appears to be primarily driven by change in the accessory gene complement, via pseudogenization, gene deletion, and horizontal gene transfer (Eppinger et al. 2006).

In a study over a more directly comparable timescale, transmission of a single clone ofEscherichia coli between the members of a family was studied for 3 years (Reeves et al. 2011). Six or more transmission events occurred, including at least two independent interspecies host jumps to the family’s dog. Amongst the 14 isolates sequenced, 20 SNPs were found, but there was no evidence of the movement of mobile elements or gene gain or loss. More of the amino acid changing mutations occurred on the lineages leading to the interspecific transmissions than expected by chance, suggesting that they may be associated with rapid adaptation to the new host.

Our data forwMelPop-CLA reflect both gene- and nucleotide-level changes. Deletion of a labile genomic region occurred rapidly, possibly beginning the process of restructuring the accessory genome. No gene gain was observed inwMelPop-CLA, as expected in a single-strain laboratory infection. Both gene loss and gene gain, however, are likely to occur inWolbachia strains over longer time periods after host jumps (Duplouy et al. 2013). The 10-bp deletion and two SNPs observed during cell line passaging all affect protein sequence and at least some are likely to have functional consequences, but there are too few changes to robustly conclude that they are due to selective processes.

The fact that we detected no further changes in the genome ofwMelPop-PGYP after 4 years in mosquitoes might indicate that the burst of substitutions in cell lines reflected adaptation to the mosquito. Alternatively, it may mean that these changes could only be fixed by drift in the permissive cell line environment. In either case, this has implications for the release ofWolbachia-infected mosquitoes for biocontrol (McGraw and O’Neill 2013): We do not expect to observe rapid evolution ofwMelPop-PGYP in released mosquitoes, meaning that pathogen-blocking and life-shortening phenotypes are unlikely to be quickly lost due to changes inWolbachia.

Use for Functional Genetics

The lack of a genetic transformation technique inWolbachia has inhibited our ability to perform functional genetics on this increasingly important bacterial species. Mutants generated during the maintenance ofWolbachia in cell culture for long periods could be isolated and exploited to create novelWolbachia infections in insects or to fine-tune existing transinfected lines using attenuated or virulent variants. Closely relatedWolbachia variants may also allow comparative genomic studies to link genotype and phenotype as an alternative to genetic transformation, for example, by comparing the phenotypes induced byWolbachia strains before and after deletion, insertion, or mutation events. The five mutations that we identified inwMelPop-CLA provide a short list of targets for further functional characterization to investigate potential mechanisms by whichWolbachia might adapt to new hosts. Furthermore, using a similar approach to understand the molecular basis forWolbachia-mediated pathogen interference could potentially open new avenues to develop novel antiviral/antimalarial compounds and to identify alternative pathways to target these pathogens.

1. Albers CA, et al. Dindel: accurate indel calls from short-read data. Genome Res. 2011;21:961–973. doi: 10.1101/gr.112326.110. [DOI] [PMC free article] [PubMed] [Google Scholar]
2. Andrews ES, et al. Reactive oxygen species production and Brugia pahangi survivorship in Aedes polynesiensis with artificial Wolbachia infection types. PLoS Pathog. 2012;8:e1003075. doi: 10.1371/journal.ppat.1003075. [DOI] [PMC free article] [PubMed] [Google Scholar]
3. Bridges CB, Brehme KS. The mutants of Drosophila melanogaster. Washington: Carnegie Institute; 1944. [Google Scholar]
4. Brownlie JC, et al. Evidence for metabolic provisioning by a common invertebrate endosymbiont, Wolbachia pipientis, during periods of nutritional stress. PLoS Pathog. 2009;5:e1000368. doi: 10.1371/journal.ppat.1000368. [DOI] [PMC free article] [PubMed] [Google Scholar]
5. Cerveau N, et al. Short- and long-term evolutionary dynamics of bacterial insertion sequences: insights from Wolbachia endosymbionts. Genome Biol Evol. 2011;3:1175–1186. doi: 10.1093/gbe/evr096. [DOI] [PMC free article] [PubMed] [Google Scholar]
6. Cook PE, McMeniman CJ, O’Neill SL. Modifying insect population age structure to control vector-borne disease. In: Aksoy S, editor. Transgenesis and the management of vector-borne disease. Austin, Texas: Landes Bioscience; 2008. pp. 126–140. [DOI] [PubMed] [Google Scholar]
7. Darling ACE, Mau B, Blattner FR, Perna NT. Mauve: multiple alignment of conserved genomic sequence with rearrangements. Genome Res. 2004;14:1394–1403. doi: 10.1101/gr.2289704. [DOI] [PMC free article] [PubMed] [Google Scholar]
8. David M, et al. SHRiMP2: sensitive yet practical short read mapping. Bioinformatics. 2011;27:1011–1012. doi: 10.1093/bioinformatics/btr046. [DOI] [PubMed] [Google Scholar]
9. Delihas N. Impact of small repeat sequences on bacterial genome evolution. Genome Biol Evol. 2011;3:959–973. doi: 10.1093/gbe/evr077. [DOI] [PMC free article] [PubMed] [Google Scholar]
10. Dobson SL, et al. Characterization of Wolbachia host cell range via the in vitro establishment of infections. Appl Environ Microbiol. 2002;68:656–660. doi: 10.1128/AEM.68.2.656-660.2002. [DOI] [PMC free article] [PubMed] [Google Scholar]
11. Duplouy A, et al. Draft genome sequence of the male-killing Wolbachia strain wBol1 reveals recent horizontal gene transfers from diverse sources. BMC Genomics. 2013;14:20. doi: 10.1186/1471-2164-14-20. [DOI] [PMC free article] [PubMed] [Google Scholar]
12. Duron O, et al. Transposable element polymorphism of Wolbachia in the mosquito Culex pipiens: evidence of genetic diversity, superinfection and recombination. Mol Ecol. 2005;14:1561–1573. doi: 10.1111/j.1365-294X.2005.02495.x. [DOI] [PubMed] [Google Scholar]
13. Ellis JC, Brown JW. The evolution of RNase P and its RNA. In: Liu F, Altman S, editors. Ribonuclease P. New York: Springer; 2010. pp. 17–40. [Google Scholar]
14. Eppinger M, et al. Who ate whom? Adaptive Helicobacter genomic changes that accompanied a host jump from early humans to large felines. PLoS Genet. 2006;2:1097–1110. doi: 10.1371/journal.pgen.0020120. [DOI] [PMC free article] [PubMed] [Google Scholar]
15. Fang G, et al. Genome-wide mapping of methylated adenine residues in pathogenic Escherichia coli using single-molecule real-time sequencing. Nat Biotechnol. 2012;30:1232–1239. doi: 10.1038/nbt.2432. [DOI] [PMC free article] [PubMed] [Google Scholar]
16. Fernandes AP, Holmgren A. Glutaredoxins: glutathione-dependent redox enzymes with functions far beyond a simple thioredoxin backup system. Antioxid Redox Signal. 2004;6:63–74. doi: 10.1089/152308604771978354. [DOI] [PubMed] [Google Scholar]
17. Garrison E, Marth G. Bayesian haplotype-based polymorphism discovery. 2012 arXiv:1207.3907. [Google Scholar]
18. Gray YHM. It takes two transposons to tango—transposable-element-mediated chromosomal rearrangements. Trends Genet. 2000;16:461–468. doi: 10.1016/s0168-9525(00)02104-1. [DOI] [PubMed] [Google Scholar]
19. Guinane CM, et al. Evolutionary genomics of Staphylococcus aureus reveals insights into the origin and molecular basis of ruminant host adaptation. Genome Biol Evol. 2010;2:454–466. doi: 10.1093/gbe/evq031. [DOI] [PMC free article] [PubMed] [Google Scholar]
20. Hannah AM. Proceedings of the 8th International Congress of Genetics (Hereditas Suppl Vol.) Stockholm: 1949. Radiation mutations involving the cut locus in Drosophila; pp. 588–589. [Google Scholar]
21. Hedges LM, Brownlie JC, O’Neill SL, Johnson KN. Wolbachia and virus protection in insects. Science. 2008;322:702–702. doi: 10.1126/science.1162418. [DOI] [PubMed] [Google Scholar]
22. Heusipp G, Falker S, Schmidt MA. DNA adenine methylation and bacterial pathogenesis. Int J Med Microbiol. 2007;297:1–7. doi: 10.1016/j.ijmm.2006.10.002. [DOI] [PubMed] [Google Scholar]
23. Hilgenboecker K, et al. How many species are infected with Wolbachia?—a statistical analysis of current data. FEMS Microbiol Lett. 2008;281:215–220. doi: 10.1111/j.1574-6968.2008.01110.x. [DOI] [PMC free article] [PubMed] [Google Scholar]
24. Hoffmann AA, et al. Successful establishment of Wolbachia in Aedes populations to suppress dengue transmission. Nature. 2011;476:454–457. doi: 10.1038/nature10356. [DOI] [PubMed] [Google Scholar]
25. Hosokawa T, et al. Wolbachia as a bacteriocyte-associated nutritional mutualist. Proc Natl Acad Sci U S A. 2010;107:769–774. doi: 10.1073/pnas.0911476107. [DOI] [PMC free article] [PubMed] [Google Scholar]
26. Ibba M, Soll D. Aminoacyl-tRNA synthesis. Annu Rev Biochem. 2000;69:617–650. doi: 10.1146/annurev.biochem.69.1.617. [DOI] [PubMed] [Google Scholar]
27. Iturbe-Ormaetxe I, Burke GR, Riegler M, O’Neill SL. Distribution, expression, and motif variability of ankyrin domain genes in Wolbachia pipientis. J Bacteriol. 2005;187:5136–5145. doi: 10.1128/JB.187.15.5136-5145.2005. [DOI] [PMC free article] [PubMed] [Google Scholar]
28. Iturbe-Ormaetxe I, et al. A simple protocol to obtain highly pure Wolbachia endosymbiont DNA for genome sequencing. J Microbiol Methods. 2011;84:134–136. doi: 10.1016/j.mimet.2010.10.019. [DOI] [PubMed] [Google Scholar]
29. Krasilnikov AS. Ribonucleoprotein ribonucleases P and MRP. In: Nicholson AW, editor. Ribonucleases. New York: Springer; 2011. pp. 319–342. [Google Scholar]
30. Leclercq S, Giraud I, Cordaux R. Remarkable abundance and evolution of mobile group II introns in Wolbachia bacterial endosymbionts. Mol Biol Evol. 2011;28:685–697. doi: 10.1093/molbev/msq238. [DOI] [PubMed] [Google Scholar]
31. Li H, Durbin R. Fast and accurate short read alignment with Burrows-Wheeler transform. Bioinformatics. 2009;25:1754–1760. doi: 10.1093/bioinformatics/btp324. [DOI] [PMC free article] [PubMed] [Google Scholar]
32. Li H, et al. The sequence alignment/map format and SAMtools. Bioinformatics. 2009;25:2078–2079. doi: 10.1093/bioinformatics/btp352. [DOI] [PMC free article] [PubMed] [Google Scholar]
33. Low DA, Weyand NJ, Mahan MJ. Roles of DNA adenine methylation in regulating bacterial gene expression and virulence. Infect Immun. 2001;69:7197–7204. doi: 10.1128/IAI.69.12.7197-7204.2001. [DOI] [PMC free article] [PubMed] [Google Scholar]
34. Lowder BV, et al. Recent human-to-poultry host jump, adaptation, and pandemic spread of Staphylococcus aureus. Proc Natl Acad Sci U S A. 2009;106:19545–19550. doi: 10.1073/pnas.0909285106. [DOI] [PMC free article] [PubMed] [Google Scholar]
35. Mackay TFC, et al. The Drosophila melanogaster genetic reference panel. Nature. 2012;482:173–178. doi: 10.1038/nature10811. [DOI] [PMC free article] [PubMed] [Google Scholar]
36. Martin F, et al. Characterization of a thermosensitive Escherichia coli aspartyl-tRNA synthetase mutant. J Bacteriol. 1997;179:3691–3696. doi: 10.1128/jb.179.11.3691-3696.1997. [DOI] [PMC free article] [PubMed] [Google Scholar]
37. McGraw EA, O’Neill SL. Beyond insecticides: new thinking on an ancient problem. Nat Rev Microbiol. 2013;11:181–193. doi: 10.1038/nrmicro2968. [DOI] [PubMed] [Google Scholar]
38. McMeniman CJ, et al. Host adaptation of a Wolbachia strain after long-term serial passage in mosquito cell lines. Appl Environ Microbiol. 2008;74:6963–6969. doi: 10.1128/AEM.01038-08. [DOI] [PMC free article] [PubMed] [Google Scholar]
39. McMeniman CJ, et al. Stable introduction of a life-shortening Wolbachia infection into the mosquito Aedes aegypti. Science. 2009;323:141–144. doi: 10.1126/science.1165326. [DOI] [PubMed] [Google Scholar]
40. Meglecz E, et al. A shot in the genome: how accurately do shotgun 454 sequences represent a genome? BMC Res Notes. 2012;5:259. doi: 10.1186/1756-0500-5-259. [DOI] [PMC free article] [PubMed] [Google Scholar]
41. Min KT, Benzer S. Wolbachia, normally a symbiont of Drosophila, can be virulent, causing degeneration and early death. Proc Natl Acad Sci U S A. 1997;94:10792–10796. doi: 10.1073/pnas.94.20.10792. [DOI] [PMC free article] [PubMed] [Google Scholar]
42. Minoche AE, Dohm JC, Himmelbauer H. Evaluation of genomic high-throughput sequencing data generated on Illumina HiSeq and Genome Analyzer systems. Genome Biol. 2011;12:R112. doi: 10.1186/gb-2011-12-11-r112. [DOI] [PMC free article] [PubMed] [Google Scholar]
43. Moreira LA, et al. A Wolbachia symbiont in Aedes aegypti limits infection with dengue, chikungunya, and Plasmodium. Cell. 2009;139:1268–1278. doi: 10.1016/j.cell.2009.11.042. [DOI] [PubMed] [Google Scholar]
44. Murray IA, et al. The methylomes of six bacteria. Nucleic Acids Res. 2012;40:11450–11462. doi: 10.1093/nar/gks891. [DOI] [PMC free article] [PubMed] [Google Scholar]
45. O’Neill SL, et al. In vitro cultivation of Wolbachia pipientis in an Aedes albopictus cell line. Insect Mol Biol. 1997;6:33–39. doi: 10.1046/j.1365-2583.1997.00157.x. [DOI] [PubMed] [Google Scholar]
46. Osborne SE, Leong YS, O’Neill SL, Johnson KN. Variation in antiviral protection mediated by different Wolbachia strains in Drosophila simulans. PLoS Pathog. 2009;5:e1000656. doi: 10.1371/journal.ppat.1000656. [DOI] [PMC free article] [PubMed] [Google Scholar]
47. Paraskevopoulos C, et al. Toward a Wolbachia multilocus sequence typing system: discrimination of Wolbachia strains present in Drosophila species. Curr Microbiol. 2006;53:388–395. doi: 10.1007/s00284-006-0054-1. [DOI] [PubMed] [Google Scholar]
48. Prevost G, et al. Study of the arrangement of the functional domains along the yeast cytoplasmic aspartyl-transfer RNA-synthetase. Eur J Biochem. 1989;180:351–358. doi: 10.1111/j.1432-1033.1989.tb14655.x. [DOI] [PubMed] [Google Scholar]
49. Rasgon JL, Styer LM, Scott TW. Wolbachia-induced mortality as a mechanism to modulate pathogen transmission by vector arthropods. J Med Entomol. 2003;40:125–132. doi: 10.1603/0022-2585-40.2.125. [DOI] [PubMed] [Google Scholar]
50. Reeves PR, et al. Rates of mutation and host transmission for an Escherichia coli clone over three years. PLoS One. 2011;6:e26907. doi: 10.1371/journal.pone.0026907. [DOI] [PMC free article] [PubMed] [Google Scholar]
51. Richardson MF, et al. Population genomics of the Wolbachia endosymbiont in Drosophila melanogaster. PLoS Genet. 2012;8:e1003129. doi: 10.1371/journal.pgen.1003129. [DOI] [PMC free article] [PubMed] [Google Scholar]
52. Riegler M, et al. Tandem repeat markers as novel diagnostic tools for high resolution fingerprinting of Wolbachia. BMC Microbiol. 2012;12:S12. doi: 10.1186/1471-2180-12-S1-S12. [DOI] [PMC free article] [PubMed] [Google Scholar]
53. Riegler M, Sidhu M, Miller WJ, O’Neill SL. Evidence for a global Wolbachia replacement in Drosophila melanogaster. Curr Biol. 2005;15:1428–1433. doi: 10.1016/j.cub.2005.06.069. [DOI] [PubMed] [Google Scholar]
54. Rouhier N, Couturier J, Johnson MK, Jacquot JP. Glutaredoxins: roles in iron homeostasis. Trends Biochem Sci. 2010;35:43–52. doi: 10.1016/j.tibs.2009.08.005. [DOI] [PMC free article] [PubMed] [Google Scholar]
55. Saridaki A, et al. Wolbachia prophage DNA adenine methyltransferase genes in different Drosophila-Wolbachia associations. PLoS One. 2011;6:e19708. doi: 10.1371/journal.pone.0019708. [DOI] [PMC free article] [PubMed] [Google Scholar]
56. Sinkins SP, Curtis CF, O’Neill SL. The potential application of inherited symbiont systems to pest control. In: O'Neill SL, Hoffmann AA, Werren JH, editors. Influential passengers. Oxford: Oxford University Press; 1997. pp. 155–175. [Google Scholar]
57. Smith HL, Rajan TV. Tetracycline inhibits development of the infective-stage larvae of filarial nematodes in vitro. Exp Parasitol. 2000;95:265–270. doi: 10.1006/expr.2000.4525. [DOI] [PubMed] [Google Scholar]
58. Stamatakis A. RAxML-VI-HPC: maximum likelihood-based phylogenetic analyses with thousands of taxa and mixed models. Bioinformatics. 2006;1:2688–2690. doi: 10.1093/bioinformatics/btl446. [DOI] [PubMed] [Google Scholar]
59. Sun LV, et al. Determination of Wolbachia genome size by pulsed-field gel electrophoresis. J Bacteriol. 2001;183:2219–2225. doi: 10.1128/JB.183.7.2219-2225.2001. [DOI] [PMC free article] [PubMed] [Google Scholar]
60. Sun LV, Riegler M, O’Neill SL. Development of a physical and genetic map of the virulent Wolbachia strain wMelPop. J Bacteriol. 2003;185:7077–7084. doi: 10.1128/JB.185.24.7077-7084.2003. [DOI] [PMC free article] [PubMed] [Google Scholar]
61. Teixeira L, Ferreira A, Ashburner M. The bacterial symbiont Wolbachia induces resistance to RNA viral infections in Drosophila melanogaster. PLoS Biol. 2008;6:2753–2763. doi: 10.1371/journal.pbio.1000002. [DOI] [PMC free article] [PubMed] [Google Scholar]
62. Treangen TJ, Salzberg SL. Repetitive DNA and next-generation sequencing: computational challenges and solutions. Nat Rev Genet. 2012;13:36–46. doi: 10.1038/nrg3117. [DOI] [PMC free article] [PubMed] [Google Scholar]
63. Valencia JI, Muller HJ. Proceedings of the 8th International Congress of Genetics (Hereditas Suppl Vol.); 1949. The mutational potentialities of some individual loci in Drosophila; pp. 681–684. [Google Scholar]
64. Walker T, et al. The wMel Wolbachia strain blocks dengue and invades caged Aedes aegypti populations. Nature. 2011;476:450–453. doi: 10.1038/nature10355. [DOI] [PubMed] [Google Scholar]
65. Werren JH, Baldo L, Clark ME. Wolbachia: master manipulators of invertebrate biology. Nat Rev Microbiol. 2008;6:741–751. doi: 10.1038/nrmicro1969. [DOI] [PubMed] [Google Scholar]
66. Woolfit M, Iturbe-Ormaetxe I, McGraw EA, O’Neill SL. An ancient horizontal gene transfer between mosquito and the endosymbiotic bacterium Wolbachia pipientis. Mol Biol Evol. 2009;26:367–374. doi: 10.1093/molbev/msn253. [DOI] [PubMed] [Google Scholar]
67. Wu M, et al. Phylogenomics of the reproductive parasite Wolbachia pipientis wMel: a streamlined genome overrun by mobile genetic elements. PLoS Biol. 2004;2:327–341. doi: 10.1371/journal.pbio.0020069. [DOI] [PMC free article] [PubMed] [Google Scholar]
68. Zug R, Hammerstein P. Still a host of hosts for Wolbachia: analysis of recent data suggests that 40% of terrestrial arthropod species are infected. PLoS One. 2012;7:e38544. doi: 10.1371/journal.pone.0038544. [DOI] [PMC free article] [PubMed] [Google Scholar]

Supplementary Materials