The expression of one ankyrin pk2 allele of the WO prophage is correlated with the Wolbachia feminizing effect in isopods

Background The maternally inherited α-Proteobacteria Wolbachia pipientis is an obligate endosymbiont of nematodes and arthropods, in which they induce a variety of reproductive alterations, including Cytoplasmic Incompatibility (CI) and feminization. The genome of the feminizing wVulC Wolbachia strain harboured by the isopod Armadillidium vulgare has been sequenced and is now at the final assembly step. It contains an unusually high number of ankyrin motif-containing genes, two of which are homologous to the phage-related pk1 and pk2 genes thought to contribute to the CI phenotype in Culex pipiens. These genes encode putative bacterial effectors mediating Wolbachia-host protein-protein interactions via their ankyrin motifs. Results To test whether these Wolbachia homologs are potentially involved in altering terrestrial isopod reproduction, we determined the distribution and expression of both pk1 and pk2 genes in the 3 Wolbachia strains that induce CI and in 5 inducing feminization of their isopod hosts. Aside from the genes being highly conserved, we found a substantial copy number variation among strains, and that is linked to prophage diversity. Transcriptional analyses revealed expression of one pk2 allele (pk2b2) only in the feminizing Wolbachia strains of isopods. Conclusions These results reveal the need to investigate the functions of Wolbachia ankyrin gene products, in particular those of Pk2, and their host targets with respect to host sex manipulation.


Background
The evolutionary success of the maternally inherited α-Proteobacteria Wolbachia pipientis is partly due to its ability to manipulate host reproduction to favour vertical transmission from mother to offspring. Wolbachia are also able to switch between hosts via horizontal transfer, which contributes to the impressive diversity and range of infected hosts [1]. These obligate endosymbionts are found in most filarial nematodes and are estimated to be present in~60% of arthropod species [2][3][4]. In arthropods, Wolbachia are considered to be sex-parasites because they alter compatibility between eggs and sperm, feminize or kill males, or induce parthenogenesis [2,5,6]. Since Wolbachia remain unculturable endosymbionts, comparative genomics and evolutionary approaches are particularly useful for identifying putative bacterial determinants involved in Wolbachia-host interactions.
Recent genome analyses of different Wolbachia strains revealed a surprisingly high number of ankyrin domaincontaining genes (ank genes) [7][8][9][10][11]. Their presence is suggested to be the result of lateral gene transfer since they are mostly found in eukaryotes but in few bacterial and viral genomes [12,13]. The 33-residue ankyrin repeats (ANK) form tandem arrays that mediate specific protein-protein interactions and have diverse functions in transcription initiation, cell cycle regulation and signalling, cytoskeleton integrity, ion transport, inflammatory responses and development [12,14]. The two closely related intracellular bacteria Anaplasma phagocytophilum and Ehrlichia chaffeensis secrete ankyrin proteins (AnkA and p200, respectively) that bind to host DNA and/or proteins [15,16]. It has been demonstrated that AnkA plays an important role in facilitating intracellular infection [17] whereas p200 is thought to affect host cell gene transcription and promote the survival of the pathogen [16]. Hence it has been suggested that ank genes encode Wolbachia effectors that alter host biology [18,19].
Several studies have suggested that Wolbachia ANK proteins were implicated in the molecular basis of Cytoplasmic Incompatibility (CI) [8,9,[20][21][22][23]. Sequence divergence between closely related Wolbachia strains causing distinct CI types in Culex pipiens mosquito populations has been found in two ank genes (pk1 and pk2), among a subset of phylogenetic markers [20,[22][23][24][25]. Thus, these polymorphic pk1 and pk2 ank genes, located within a socalled WO prophage region of Wolbachia genome, are suggested to contribute to the CI phenotype. Consistent with this argument, expression of the pk2 gene occurred specifically in female mosquitoes [8,22,23]. Moreover, a premature stop codon was found in the pk2 gene of the Wolbachia strain (wAu) that is unable to cause CI in D. simulans [21].
In this study, we aimed to determine whether the prophage pk1 and pk2 ankyrin genes were involved in the CI phenotype described in three Wolbachia-infected species of terrestrial isopods. We also investigated whether these genes were conserved and expressed in Wolbachia strains inducing feminization, the main Wolbachia phenotype described for this group of hosts [2]. From the genome of the feminizing wVulC Wolbachia strain that infects the isopod Armadillidium vulgare (the genome completion is currently being done by our group in the frame of the European Wolbachia project: EuWol), we annotated the pk1 and pk2 alleles among all ank genes identified from the wVulC contigs. We investigated the distribution, copy number and expression patterns of both genes in seven additional Wolbachia strains that induce either CI or feminization in isopods. We identified a large copy number variation of the pk1 and pk2 genes among Wolbachia strains, which is probably coupled to prophage evolution. Surprisingly, our results also revealed that expression of one pk2 allele (pk2b2) is only detected in feminizing Wolbachia strains and never in the three CI-inducing strains of isopods.

Results
Characterization and distribution of pk1 and pk2 genes Six copies of the pk1 gene and three copies of the pk2 gene were identified in the contig assembly of the wVulC genome ( Table 1). Each of the six putative prophage regions of the assembly contains one pk1 allele and three of these prophages also harbour one pk2 allele (Table 1). Two wVulC pk1 alleles (ANK46a/b and ANK60a/b) and one pk2 allele (ANK40a/b) were each found in two identical copies. These results were confirmed by Southern blotting (Additional file 1: Figure S1) and are consistent with the sequencing of PCR products (Table 1).
Relationships between these alleles and all the published sequences of Wolbachia pk1 and pk2 genes were assessed based on the DNA sequences encoding ANKrepeats. No evidence of recombination events has been detected in the alignments of pk1 and pk2 gene regions encoding ANK motifs (see Methods). In the pk1 sequences, the number of variable sites is 467 out of 1068, of which 408 are informative. Similarly, there are 66 informative sites in the 292 bp-long pk2 sequence alignment. The resulting genetic networks show that pk1 and pk2 sequences group in two clusters ( Figure 1). Based on these relationships, we defined two different types of pk1 genes (named pk1a and pk1b) and two different types of pk2 genes (named pk2a and pk2b) ( Figure 1). The wVulC genome contains a single copy of the pk1a type (ANK25) and 5 copies belonging to the pk1b type: ANK32, ANK46a/b and ANK60a/b ( Figure 1A). The wVulC pk1a gene clusters together with the pk1a type gene of wMel ( Figure 1A, 15.9% divergent) whereas it shares 55.2 to 65.2 % identity with the five wVulC pk1b type sequences. The region encoding ANK repeats in the wVulC ANK60a/b alleles is related to the pk1 sequence from the wCauB2 prophage ( Figure 1A). wVulC ANK46a/b alleles are closely related to the pk1 gene from the wCauB3 prophage ( Figure 1A). ANK32 seems somehow related to the pk1 gene from the WOVitA1 prophage ( Figure 1A). The wVulC genome also harbours three genes of the pk2b type, ANK48 and ANK40a/b, further called pk2b1 and pk2b2 alleles respectively. In contrast to the pk1 gene, all three wVulC pk2 alleles form a cluster in the gene network ( Figure 1B). Their closest relative is the pk2 gene harboured by the WOVitA4 prophage of the Wolbachia strain endosymbiont of Nasonia vitripennis ( Figure 1B).
Using the same primer set as for wVulC (Additional file 1: Table S1), the taxonomic distribution of pk1 and pk2 genes was extended by PCR to seven Wolbachia strains that induce either CI or feminization in isopods. All these strains of isopods are known to belong to the B-supergroup of Wolbachia whatever the phylogenetic marker used [2]. They do not form separate monophyletic clades according to the phenotype they induce in their hosts based on the wsp gene (Additional file 1: Figure S2). We also investigated the copy number variation by Southern blot analyses of EcoRI or BamHI digested DNA using pk1a, pk1b and pk2b1 probes which, according to sequence identities, preferentially hybridized on pk1a, pk1b and pk2b types, respectively (Table 2 & Additional file 1: Figure S1). In congruence with amplification and sequencing data, the pk1a and pk1b probes revealed two to six copies of the pk1 gene in the studied strains ( Table 2). By direct sequencing of the PCR products, we found that the pk1a gene of Wolbachia strains of C. convexus, P. pruinosus, A. vulgare (wVulM) and A. nasatum harboured 1, 1, 2 and 3 EcoRI sites, respectively, explaining the discrepancy between the number of bands observed by Southern blots, and the number of different sequences obtained (Table 2 & Additional file 1: Figure S1). Similarly, two pk1b alleles of the Wolbachia strain of A. nasatum contained one BamHI restriction site. Each of the two more intense Southern Blot signals (Additional file 1: Figure S1) revealed the presence of two identical copies wVulC pk1b alleles, as confirmed by the analysis of contigs. Furthermore, Southern blots using a pk2b1 probe in combination with sequencing data revealed three copies of the pk2 gene in all strains tested except one ( Table 2 & Additional file 1: Figure S1). In the Wolbachia strain of P. pruinosus, sequences of PCR products revealed two identical pk2 alleles, each containing one BamHI restriction site explaining the five signals obtained by Southern blotting ( Table 2 & Additional file 1: Figure S1). Moreover, no signal was obtained from digested and undigested DNA of Wolbachia-free ovaries of isopod (non-infected population from Nice, France), which confirmed the Wolbachia origin of the pk1 and pk2 genes.
terminal region but a variable number of ANK motifs ranging from 8 to 10 (Additional file 1: Figure S3). In wVulC, ANK46a/b and ANK60a/b sequences (pk1b type) are shorter in their N-terminal region than the other Pk1 translated sequences (42 and 62 amino acids, respectively). One indel at position 117 of the DNA sequence of wVulC ANK46a/b is responsible for a frame shift, which splits the gene into two ORFs homologous to the fulllength pk1 of other strains. ANK60a/b sequences are shortened by a transposase gene insertion in the 5′ region. In contrast, pk2 translated sequences are more conserved (84.5 to 100% identity) among Wolbachia strains than pk1. All Pk2 amino acid sequences harbour 3 ANK motifs except in the wAu strain (host: D. simulans) in which a premature stop codon disrupts the third motif (Additional file 1: Figure S3).
Comparative analysis of pk1 and pk2 mRNA expression in CI and feminizing Wolbachia strains RT-PCR using allele-specific primers was performed to examine the expression patterns of pk1 and pk2 mRNA in adult gonads of isopods harbouring CI-inducing or feminizing Wolbachia strains (Figure 2). Evidence of expression was observed for all copies of pk1 and pk2 genes except for one allele of the pk2b type (Figure 2A). Although the pk2 alleles were amplified by PCR using DNA of all the three isopod CI-inducing Wolbachia strains (positive controls only shown for wVulC), RT-PCR did not yield any detectable products in these strains using the pk2b2 allele-specific primers ( Figure 2B). In contrast, the pk2b2 allele was clearly expressed in all the feminizing Wolbachia strains ( Figure 2B). In hosts where both males and females are infected by CI-inducing or feminizing strains, no clear sex-specific differences were observed in pk1 and pk2 expression (Figure 2A). We further examined the expression of pk2b2 and another prophage gene, orf7 which encodes the phage capsid, in several tissues of A. vulgare females harbouring the feminizing wVulC strain ( Figure 2C). While orf7 was expressed only in ovaries, the host tissue where the density of Wolbachia is higher, transcription of pk2b2 was revealed in all tissues tested (except the brain) ( Figure 2C).

Discussion
In this study, we found that a large copy number variation of pk1 and pk2 genes exists among Wolbachia strains, which is probably coupled to prophage dynamics and evolution. Copy number divergence in the ankyrin pk1 and pk2 is consistent with the results of previous Southern blotting analyses using the minor capsid orf7 phage gene [28]. Four different orf7 paralogs had already been identified in the wVulC strain through cloning and sequencing of heterogeneous PCR products [28]. Since multiple infections of Wolbachia in a single individual have never been observed in isopods, we can conclude that the phage WO is likely to be present in several copies in each Wolbachia strain. Our observations of Wolbachia strains of isopods suggest that dynamics of the prophage pk1 and pk2 genes is similar to that observed in the wRi and wPip-Pel genomes [8,9]. While the pk1 and pk2 ank genes are each found as a single copy within the WO-B prophage of the wMel genome [11], sequencing of the CI-inducing wPip-Pel, wRi and wVitA Wolbachia prophages revealed a duplication of these genes together with the WO-B-like prophages [8,9,27]. For instance, in the wPip-Pel genome, the three pk1 and the three pk2 genes are spread among the five different prophages which are closely related to the WO-B wMel prophage [8]. Hence, the divergence in the pk1 and pk2 gene copy number between Wolbachia strains may be explained by mechanisms related to bacterial genome organization and modulation of gene copy number [26,[29][30][31][32]. As an example, two pseudogenes (wRi_ANK29 and ANK31) out the four copies of the pk1 gene in wRi, are spread in the WORiB prophage (previously annotated WO-C prophage [9], see Table 1) and may have originally been a single pk1 gene further disrupted by an insertion sequence ISWpi7. On the other hand, the high GC content of pk2 supports the occurrence of recent lateral transfers of prophage fragments containing the pk2 gene but not necessarily pk1 in the Wolbachia genomes. However, we cannot exclude the hypothesis that linkage disequilibrium occurs between pk1 and pk2 genes that are separated by at least 6.7 kilobases, representing less than 0.04% of the whole genome size. These results also highlight the genomic plasticity The number of pseudogenes is indicated in brackets. ND: non determined. of the prophage region among Wolbachia strains as part of the global plasticity observed in the Wolbachia genomes [33]. Maintenance of such "mobile elements" in Wolbachia strains of arthropods may be due to the absence of, or a reduced efficiency of selection on the prophages. Nevertheless, the purifying selection acting on these pk1 and pk2 genes suggest that maintenance of sequences confers an adaptive advantage. Besides identifying mosaic prophages, our results also reveal the differential expression of one pk2 ankyrin according to the Wolbachia phenotype they induce (CI vs. feminization). One allele (pk2b2) is only expressed in the feminizing strains and never in the three CI-inducing strains of isopods. In contrast to the observations for wPip [22,23], expression pattern of pk2b2 suggests that this allele is not involved in CI in isopods. In two recent studies, it has been shown that expression of pk1 and pk2 genes from wMel was not correlated with the CI phenotype in D. melanogaster [34,35]. Our transcriptional result rather leads to the hypothesis that this pk2b2 allele is involved in the feminization of isopod hosts. This hypothesis is strengthened by the observation that the pk2b2 allele is expressed in all A. vulgare tissues (except in the brain) whereas another prophage gene (orf7) is only expressed in ovaries. Furthermore, no differential expression of pk1 and pk2 genes was identified between sexes in isopods when either CI-inducing or feminizing Wolbachia infects both males and females. This result differs from those of Sinkins and colleagues who showed that in some CI-inducing wPip variants, the three pk2 genes (the two identical wPip_ANK12 and wPip_ANK25, and wPip_ANK16) are highly expressed in females but never in males [22,23].
Our data do not enable us to explain why pk2b2 is only expressed in feminizing strains of Wolbachia whereas its homologs, also found in CI-inducing strains, are associated with CI phenotype in mosquitoes [22,23]. First, one can suggest that this allele has been inactivated or importantly down regulated in the CI-inducing strains of isopods. Change in regulatory element repertoire and divergence in patterns of expression may occur after small-scale duplication of the genome [36]. A corollary to a change in location, paralogous and homologous pk2 copies within and among Wolbachia strains would have followed different evolutionary trajectories leading to such a phenotypic diversity. Second, genomic imprinting, process by which genes are expressed from only one parental allele due to epigenetic mechanism, can be considered as a molecular mechanism underlying the diversity of phenotypes. Recently, early changes in gene imprinting and aberrant expression of specific genes have been shown to be coupled to parthenogenesis in mice embryos [37]. Third, one can suggest that genes in the pk2 family could have diverse functions. In this way, post-transcriptional modifications and dosage of Wolbachia products, as well as genetic control by the host, cannot be dismissed. As previously suggested [38], differences in Wolbachiainduced feminization as well as the presence of the bacteria in O. asellus males, may simply result from differences in bacterial dosage or in host targets. The basic molecular mechanisms that mediate Wolbachia feminization are also still unknown although it is unlikely that this effect is driven by only one gene. In A. vulgare, Wolbachia effectors may target the proteinaceous androgenic hormone or its receptor, or another major sex determinant, thereby inhibiting the androgenic gland differentiation and preventing the androgenic hormone from reaching the target tissues such as gonads and tegumental epithelium [2,39,40]. This hypothesis suggests a late action of feminizing Wolbachia on host target(s) during its development, as opposed to the very early action of other Wolbachia strains that induce parthenogenesis, CI or male killing [5,41].

Conclusions
Our results highlight a large copy number variation of both pk1 and pk2 genes among strains, likely linked to prophage diversity, and also the specific expression of one pk2 allele only in the feminizing Wolbachia strains of isopods. This correlation supports the hypothesis that phenotype-related effectors or specific strain determinants in Wolbachia are likely to be encoded by prophage genes, ankyrin-repeat encoding genes, and predicted genes of unknown function [42]. Our results thus reveal the need to search for host molecules targeted by Wolbachia ankyrins and their functions with respect to host sex manipulation by Wolbachia.
Total DNA was extracted from male and female gonads of all isopod species as described previously [48]. Infection status of each individual was confirmed by a PCR-assay based on the bacterial 16S rDNA gene using Wolbachia-specific primers (Additional file 1: Table S1) [49].

Distribution of pk1 and pk2 genes
The genome of the feminizing wVulC Wolbachia strain is at the final assembly step (whole-genome shotgunsequencing project: European Wolbachia EuWol (contract QLK3-CT2000-01079, coordinated by K. Bourtzis, University of Ioannina, Greece). This includes phage contigs of which sequences are homologous to the Wolbachia WO prophage. Annotation of the pk1 and pk2 genes was performed by protein and DNA homology searches with BLASTP and BLASTN programs [50] using the wPip-Pel pk1 and pk2 alleles as queries (see Table 1). Ankyrin and other functional motif predictions were performed by the SMART web server [51] on protein sequences.
Specific primers were designed to amplify full-length or 200-500 bp fragments of the wVulC pk1 and pk2 alleles using a standard PCR protocol as previously described (Additional file 1: Table S1) [52]. The purified PCR products were directly sequenced on both strands on an ABI PRISM 3100 Genetic Analyzer using Big Dye Terminator v3.1 Cycle Sequencing Kit (Applied Biosystems) according to the manufacturer's instructions.
pk1 or pk2 copy number variation among Wolbachia strains was assessed by Southern blotting. About 15 μg of total DNA were digested at 37°overnight with EcoRI or BamHI enzymes that did not cut any of the wVulC pk1 and pk2 alleles. Digested DNA as well as undigested DNA from non-infected ovaries used as controls (data not shown) was electrophoresed on 0.8% agarose gels and blotted to nylon membranes. Probes were obtained by PCR amplification of the wVulC full-length pk1 (pk1a and pk1b types) and pk2 (pk2b type) ank genes (Additional file 1: Table S1), labelled using [α-32 P]-dCTP by the random primer method and hybridized overnight to membranes. The final wash was performed at 52°in 0.1X SSC. Hybridized blots were imaged and analyzed using a PhosphoImager (Molecular Dynamics, Sunnyval, CA, USA).

Sequence analyses of pk1 and pk2 genes
Homologous sequences of both genes were first aligned in the server-based program MAFFT (http://align.bmr. kyushu-u.ac.jp/mafft/online/server/) using automatic settings. The resulting matrices were manually checked according to the predicted amino acid translation using BioEdit v7.0.9 [53]. MEGA5 software [54] was used to calculate nucleotide sequence divergence. For each locus, the GC content, the number of variable sites and the level of nucleotide diversity per site (Pi) were calculated. Ka/Ks likelihood analysis was also performed using the Selecton web server [55].
Recombination analysis was performed with RDP version v3.42 [56] using an alignment of non-redundant pk1 and pk2 nucleotide region encoding ANK-repeat domains. The parameters were set as follows: sequences were considered linear, the highest acceptable P value cut-off was 0.01, a Bonferroni correction was applied, consensus daughter sequences were found, gaps were included, different window sizes of variable sites were tested and 1,000 permutations were performed.
The best-fitted model of DNA evolution was estimated with jModelTest v0.1.1 [57] according to the corrected Akaike Information Criterion [58]. The selected model was TIM + G for pk1 and HKY + I for the pk2 locus encoding the ANK domain cluster. Gene genealogies were constructed using MrBayes v3.1.2 software [59,60] and supported by Bayesian and Maximum likelihood (ML) probabilities. Two Metropolis-coupled Markov chain Monte Carlo (MCMC) analyses were run for 5,000,000 generations and sampled every 250 generations. The first 25% of sampled trees were considered burn-in trees and were discarded before constructing a 50% majority rule consensus tree. ML analyses were carried out in PhyML 3.0 [61]. Node support came from 1,000 multiparametric bootstrap replicates. The networks were visualized with FigTree v1.3.1 (http://tree.bio.ed.ac.uk/software/figtree). The network tree of the wsp gene was built following an identical Bayesian methodology (model: TPM3uf + I + G) (Additional file 1: Figure S2).

Expression of ankyrin genes
Total RNAs were isolated from 20 to 50 gonads dissected from all species using the RNeasy Mini Kit (Qiagen) according to the manufacturer's instructions. Ovaries were used in A. vulgare and A. nasatum where only females are infected. After a treatment with DNaseI (2U/μL, Ambion) at 37°for 30 min, 1 μg of RNA was used for reverse transcription using Superscript III kit (Invitrogen) as described by the manufacturer. To determine the expression of each gene, 1 μL of the reverse transcriptase reaction was used as template for the RT-PCR experiments. Control of the RT reactions was performed by omitting reverse transcriptase in the negative (RT-) controls and by testing the expression of the Wolbachia 16S rDNA gene (Additional file 1: Table S1). Genomic DNA of all species was also used as a positive control of the PCR reactions as well as the one of the uninfected population (Nice, France) as negative control. Transcriptional analyses of pk2b2 and orf7 genes in several tissues of A. vulgare harbouring the feminizing wVulC Wolbachia strain were run as previously described [52]. The 16S rDNA gene expression was again used to check the quality of the extracted RNA in all these tissues.