Variation in Wolbachia cidB gene, but not cidA, is associated with cytoplasmic incompatibility mod phenotype diversity in Culex pipiens.

Endosymbiotic Wolbachia bacteria are, to date, considered the most widespread symbionts in arthropods and are the cornerstone of major biological control strategies. Such a high prevalence is based on the ability of Wolbachia to manipulate their hosts' reproduction. One manipulation called cytoplasmic incompatibility (CI) is based on the death of the embryos generated by crosses between infected males and uninfected females or between individuals infected with incompatible Wolbachia strains. CI can be seen as a modification-rescue system (or mod-resc) in which paternal Wolbachia produce mod factors, inducing embryonic defects, unless the maternal Wolbachia produce compatible resc factors. Transgenic experiments in Drosophila melanogaster and Saccharomyces cerevisiae converged towards a model where the cidB Wolbachia gene is involved in the mod function while cidA is involved in the resc function. However, as cidA expression in Drosophila males was required to observe CI, it has been proposed that cidA could be involved in both resc and mod functions. A recent correlative study in natural Culex pipiens mosquito populations has revealed an association between specific cidA and cidB variations and changes in mod phenotype, also suggesting a role for both these genes in mod diversity. Here, by studying cidA and cidB genomic repertoires of individuals from newly sampled natural C. pipiens populations harbouring wPipIV strains from North Italy, we reinforce the link between cidB variation and mod phenotype variation fostering the involvement of cidB in the mod phenotype diversity. However, no association between any cidA variants or combination of cidA variants and mod phenotype variation was observed. Taken together our results in natural C. pipiens populations do not support the involvement of cidA in mod phenotype variation.

INTRODUCTION

Wolbachia are maternally inherited endosymbiotic bacteria commonly found in arthropods and filarial nematodes (Ferri et al., 2011; Taylor, Bandi, & Hoerauf, 2005; Werren, Baldo, & Clark, 2008). More than 40% of terrestrial arthropod species are thought to be infected (Weinert, Araujo‐Jnr, Ahmed, & Welch, 2015; Zug & Hammerstein, 2012). The pervasiveness of this bacterial genus is mostly attributed to its ability to manipulate host reproduction, facilitating its spread within arthropod populations (Rousset & Raymond, 1991; Turelli & Hoffmann, 1991; Werren et al., 2008). The most commonly described Wolbachia‐induced phenotype in arthropods is cytoplasmic incompatibility (CI; Werren et al., 2008). CI is a form of conditional sterility resulting in embryonic lethality when infected males mate with uninfected females or with females infected with a different, incompatible Wolbachia strain (Atyame, Duron, et al., 2011; Bonneau, Landmann, et al., 2018; Bordenstein, O'Hara, & Werren, 2001; Breeuwer & Werren, 1990; Callaini, Riparbelli, Giordano, & Dallai, 1996; Duron et al., 2006; Laven, 1967; O'Neill & Karr, 1990). In C. pipiens where all males and females are infected, CI may be unidirectional (crossing is compatible in one direction but incompatible in the other) or bidirectional (crosses in both directions are incompatible; Atyame et al., 2014; Dumas et al., 2013; Laven, 1967; Rasgon & Scott, 2003; Sicard, Bonneau, & Weill, 2019). CI can be seen as a toxin‐antidote model or modification‐rescue model (mod‐resc) in which the Wolbachia present in the male produce a toxin (mod factors) during spermatogenesis which induces CI through embryonic defects after fertilization unless the Wolbachia present in the eggs produce compatible antidotes (resc factors) (Hurst, 1991; Poinsot, Charlat, & Merçot, 2003; Werren, 1997). Both sterile insect and pathogen blocking Wolbachia‐based methods to fight against arthropod pests and vectors rely on the ability of Wolbachia to induce CI. Knowledge on CI diversity in mosquito is required to find the better Wolbachia‐mosquito associations to optimize the success of biological control (Flores & O'Neill, 2018; Sicard et al., 2019).

Recent works have implicated the syntenic Wolbachia genes cidA and cidB, from wMel and wPip‐Buckeye strains infecting Drosophila melanogaster and Culex pipiens respectively, in CI (Figure 1; Beckmann, Ronau, & Hochstrasser, 2017; LePage et al., 2017). These genes exhibit many typical features of toxin‐antidote models (Beckmann et al., 2019a, 2019b, 2017). In both transgenic yeasts and flies, cidA Pip and cidB Pip genes were proposed to encode interacting proteins acting in a toxin‐antidote fashion with the toxicity of CidB Pip being rescued by the expression of CidA Pip (Beckmann et al., 2017). However, both cidA Mel and cidB Mel were required to induce CI in transgenic D. melanogaster (LePage et al., 2017; Shropshire & Bordenstein, 2019) and the expression of cidA Mel in transgenic D. melanogaster females was necessary and sufficient to resume the resc function (Figure 1; Shropshire & Bordenstein, 2019; Shropshire, On, Layton, Zhou, & Bordenstein, 2018). Based on these findings a two‐by‐one model was proposed, in which cidA Mel acts as a mod factor when expressed in males and as a resc factor when expressed in females, while cidB Mel acts only as a mod factor in D. melanogaster (Shropshire & Bordenstein, 2019; Shropshire et al., 2018).

Figure 1

Summary of transgenic experiments and natural population studies conducted on cidA and cidB. (a) Hypothetical representation of the portion of prophage WO containing cidA and cidB CI genes in wMel and wPip Wolbachia strains. The genome of wMel contain only one copy of cidA and cidB genes (Lindsey et al., 2018). Genomes of all the wPip strains from several wPip groups investigated contain several different copies of cidA and cidB genes (Bonneau, Atyame, et al., 2018). The deubiquitinase (DUB) region is the catalytic domain of CidB protein (Beckmann et al., 2017). This region is conserved between the cidB Pip variants (Bonneau, Atyame, et al., 2018). (b) transgenic expression of cidA and cidB from either wMel or wPip‐Buckeye genomes in Drosophila melanogaster flies. The expression of both cidA Mel and cidB Mel in D. melanogaster males is required to induce CI, neither cidA Mel or cidB Mel alone can induce CI (LePage et al., 2017; Shropshire & Bordenstein, 2019). The expression of cidA Mel alone in D. melanogaster females is sufficient to rescue CI (Shropshire et al., 2018). A two‐by‐one model of CI was proposed in D. melanogaster in which cidA Mel acts as a mod factor when expressed in males and as a resc factor when expressed in females (Shropshire & Bordenstein, 2019; Shropshire et al., 2018). The production of viable transgenic male flies expressing only cidB Pip was not possible suggesting a toxic effect of the CidB Pip protein (represented by a skull) while flies expressing both cidA Pip and cidB Pip were viable and capable of CI induction. The male flies expressing cidA Pip and cidB Pip with a disrupted catalytic DUB domain were not capable to induce CI suggesting that the DUB region is functionally involved in CI. The current transgenic data in Culex pipiens support a toxin‐antidote model where cidB Pip would encode a toxin involved in the mod function while cidA Pip would encode the antidote involved in the resc but also prevent the producer from the toxicity of CidB Pip protein. (c) cidA Pip and cidB Pip variants repertoires in natural populations of C. pipiens infected with wPipIV strains. Full variant names are not shown (they all belong to the group wPipIV) and only the letter/number of the variant appear. For instance cidA α referring to cidA_IV(α) and cidB a2 to cidB_IV(a/2). The number of pairs of genes as well as their disposition in the genome might not reflect the reality as these informations are still under investigation. Males from North African and North Italian natural populations are either compatible or incompatible with females from the Tunis isofemale line depending on the wPipIV strain they carry. All the wPipIV strains carry several cidA and cidB variants inside their genomes. The variants cidB_IV(2) [i.e., cidB_IV(a/2) and cidB_IV(b/2)] were found associated with the incompatible mod phenotype in both geographical areas while the variant cidA_IV(δ) was found associated with the incompatible mod phenotype only in the North African population. Furthermore, no other cidA variant or combination of variants was found associated with mod phenotype variation in North Italian populations suggesting that only cidB plays a role in the mod phenotype variations in C. pipiens. Finally, the cidA_IV(α) variant was detected in all the wPipIV strains regardless of their mod phenotype and their geographical origins. The ubiquity of the cidA_IV(α) variant could be responsible for the reciprocal compatibility always observed between mosquitoes infected with different wPipIV strains and suggests a role of cidA in the resc function [Colour figure can be viewed at http://www.wileyonlinelibrary.com]

Cytoplasmic incompatibility in Culex pipiens mosquitoes is characterized by its unprecedented diversity of compatibility and incompatibility relationships that is based on the diversity of the Wolbachia strains infecting this species (Atyame et al., 2014; Duron et al., 2006; Laven, 1967). All the Wolbachia infecting C. pipiens belong to the monophyletic group wPip inside the supergroup B and are closely related to wBol, the Wolbachia strain infecting the butterfly Hypolimnas bolina (Atyame, Delsuc, Pasteur, Weill, & Duron, 2011; Bleidorn & Gerth, 2018). The wPip group is divided into five groups wPipI‐V and mosquitoes infected with Wolbachia from two different wPip groups are more likely to be incompatible than mosquitoes infected with wPip strains from the same wPip group (Atyame et al., 2014). An analysis conducted on multiple crosses showed that each wPip genome must contain several mod and resc factors to account for the diversity of CI phenotypes in C. pipiens (Atyame, Duron, et al., 2011; Nor et al., 2013). These multiple mod and resc factors could theoretically be encoded by different copies (i.e., variants) of the same mod/resc genes or different mod and resc genes within the same wPip genome as already proposed in some CI models (Atyame, Duron, et al., 2011; Nor et al., 2013; Poinsot et al., 2003). To investigate the genetic basis of the unprecedented diversity of CI phenotypes found in C. pipiens mosquitoes we studied the cidA and cidB genes of wPip Wolbachia strains belonging to the groups wPipI, II, III and IV (Bonneau, Atyame, et al., 2018). All these wPip genomes exhibited several polymorphic copies of cidA Pip and cidB Pip genes (Bonneau, Atyame, et al., 2018; Bonneau, Landmann, et al., 2018). To be responsible for the CI phenotype diversity observed in C. pipiens, cidA Pip and cidB Pip genes must have different sequences in wPip strains that induce different CI phenotypes (i.e., showing different incompatibility relationships when crossed with individuals harbouring other strains). We analyzed the cidA Pip and cidB Pip variant repertoires of wPip strains showing different CI phenotypes and showed that they exhibited distinct cidA Pip and cidB Pip variants in their genomes supporting the implication of these two genes in CI phenotype diversity in C. pipiens (Bonneau, Atyame, et al., 2018).

The putative roles of the cidA Pip and cidB Pip genes in CI phenotype diversity have been further investigated by analyzing variation of cidA Pip and cidB Pip in isofemale lines infected with Wolbachia from the wPipIV phylogenetic group and exhibiting well differentiated mod phenotypes (Atyame et al., 2015; Bonneau, Atyame, et al., 2018). In C. pipiens, difference in mod phenotypes refers to differing compatibility between crosses involving focal studied males (i.e., from which we aim to infer the mod phenotype) and females infected with reference wPip strains (Bonneau, Atyame, et al., 2018). As such crosses only differ from each other because of the wPip strain harboured by the males, the difference in compatibility (i.e., egg‐raft hatching vs. no egg‐raft hatching) between these crosses allows us to qualify the mod phenotype. This way, we demonstrated that mosquitoes infected with wPipIV present two different mod phenotypes: the incompatible mod phenotype, in which males cannot produce viable progenies with females infected with wPip from group I, II or III, and the compatible mod phenotype, in which males produce viable progenies with such females (Atyame et al., 2014, 2015; Bonneau, Atyame, et al., 2018). The association of a specific cidB Pip variant (variant cidB_IV[a/2]) with the incompatible mod phenotype in 180 isofemale lines (Figure 1), together with the results of functional studies in D. melanogaster (Beckmann et al., 2017; LePage et al., 2017), demonstrated the involvement of cidB Pip in mod function and in the diversity of the CI mod phenotype in C. pipiens (Bonneau, Atyame, et al., 2018). However, the implication of cidA Pip in CI mod and/or resc function remained unclear in C. pipiens. Indeed, an association was found between a specific cidA Pip variant (cidA_IV[δ]) and the incompatible mod phenotype, suggesting a possible role for cidA Pip in mod phenotype diversity (Figure 1) (Bonneau, Atyame, et al., 2018). However, the ubiquitous presence of cidA_IV(α) in all the wPipIV strains might account for reciprocal compatibility (systematic rescue) between them suggesting a putative role for cidA Pip in the resc function (Figure 1) (Bonneau, Atyame, et al., 2018). Overall, the diversity of CI phenotypes described in C. pipiens might result from differential interactions between specific CidA and CidB variants. We have already proposed a zone of interaction between CidA and CidB proteins located in the polymorphic regions shown to be correlated with CI variation (Bonneau, Atyame, et al., 2018). It is thus possible that the association of the cidA_IV(δ) variant with the wPipIV incompatible mod phenotype results from codiversification with the cidB_IV(a/2) variant for binding adjustment, as predicted for toxin‐antidote systems.

Here, we further investigated the link between cidA and CI phenotype diversity, by studying the cidA Pip/cidB Pip variant repertoires in wPipIV strains from newly‐sampled natural populations with possibly different evolutionary histories. We screened the mod phenotype variation in four North Italian populations presenting a mixture of wPip group IV strains inducing either compatible or incompatible mod phenotype. The cidA_IV(α), thought to be associated with intra‐wPipIV compatibility was present in all the lines, regardless of their mod phenotype suggesting the putative implication of cidA Pip in the resc (antidote) function. The presence of the cidB_IV(a/2) variant in all incompatible isofemale lines supports the role of cidB in mod phenotype variation. In the analyzed North African natural populations approximately 5% of isofemale lines were qualified as incongruent as they displayed the compatible mod phenotype but carried the cidA Pip and cidB Pip variants associated (i.e., in the 95% other cases) with incompatible mod phenotypes (Bonneau, Atyame, et al., 2018). We could not study these incongruent lines previously as they were not alive anymore when we studied their cidA and cidB gene repertoires (Bonneau, Atyame, et al., 2018). In this study, we managed to study one incongruent line (i.e., being compatible while harbouring cidB_IV[a/2] variant). This line exhibited both a unique cidB repertoire and a lower expression of the cidB variant associated with incompatible mod phenotype that could contribute to explain such incongruence. Most importantly, no specific cidA Pip variant or combination of cidA Pip variants was associated with either incompatible or compatible mod phenotypes, pointing towards the absence of consequence of cidA Pip variation on CI mod diversity in C. pipiens. Overall, our data do not reject the two‐by‐one model of CI but have nothing to support it in C. pipiens. However, the variations of cidB that match changes in mod and the ubiquity of a cidA variant between compatible strains fit the expectation of a classic toxin‐antidote model.

MATERIALS AND METHODS

Mosquito collection and the construction of isofemale lines

Culex pipiens larvae and pupae were collected from four natural breeding sites in North Italy in 2017 (Roveré della Luna, San Michele all'Adige, Zambana and Mezzocorona sites, Data S1) and reared to adulthood in the laboratory. Females were then fed turkey blood (bcl Wholly Wild World) with a Hemotek membrane feeding system (Discovery Workshops, UK), to enable them to lay eggs, from which isofemale lines were established. Each egg‐raft (containing 100–250 eggs) was individually isolated for hatching, and the Wolbachia group present was determined by performing the pk1 PCR‐RFLP test (Altinli, Gunay, Alten, Weill, & Sicard, 2018) on two first‐instar larvae (L1) after extracting DNA using an acetyl trimethylammonium bromide (CTAB) protocol (Rogers & Bendich, 1989). Isofemale lines were created by rearing the offspring resulting from a single egg‐raft (thus from a single female). We established 67 isofemale lines for this study (Data S1). Isofemale lines were reared in 65 dm3 screened cages in a single room maintained at 22–25°C, under a 12 hr light/12 hr dark cycle. Larvae were fed with a mixture of shrimp powder and rabbit pellets, and adults were fed on honey solution.

Determination of CI phenotypes

CI phenotype of the isofemale lines resulting from field collection in North Italy

To be able to associate cidA and cidB variants with CI phenotype variation, the CI mod phenotype of each of the 67 North Italian isofemale lines was determined. CI mod phenotypes were characterized by crossing in the same cage males (25–50 virgin males) from each of the studied isofemale lines with females (25–50 virgin females) from the Tunis laboratory isofemale line infected with a Wolbachia strain from the wPipI group (Table S1; Duron et al., 2005). After five days in cages, the females were fed a blood meal and, five days later, egg‐rafts were collected and deposited into 24‐well plates. The CI mod status of each cross was determined by assessing eggs‐raft hatching status. All unhatched eggs‐rafts were checked for fertilization by observing embryonic development with a light microscope (Axiophot2, Zeiss), as described by Duron and Weill (2006). Two type of crosses were found: crosses with only fertilized unhatched eggs‐rafts which were qualified as incompatible and crosses with only hatched eggs‐rafts which were qualified as compatible. No crosses resulting in both fertilized hatched and fertilized unhatched eggs‐rafts were found. Thus isofemale lines in which the males were incompatible with females from the Tunis line were described as incompatible isofemale lines, whereas isofemale lines in which the males were compatible with females from the Tunis line were described as compatible isofemale lines.

Capacity of Michele26 line to induce CI

Because Michele26 males were not able to induce CI when crossed with Tunis females, the capacity of Michele26 to induce CI at all was tested. The capacity of Michele26 line to induce CI was tested by crossing 25 virgin males with 50 virgin uninfected females from the SlabTC laboratory line. SlabTC line was obtained from the Slab laboratory line treated with tetracycline as described in Duron et al. (2006).

Reciprocal compatibility of isofemale lines infected with wPipIV strains

To have a better support of the hypothesis that the cidA_IV(α) variant might be associated with the reciprocal compatibility of the isofemale lines infected with wPipIV, crosses between males and females infected with different wPipIV strains all harbouring cidA_IV(α) variant were performed. The reciprocal compatibility from seven laboratory wPipIV infected lines showing different cidA‐cidB repertoires but all exhibiting cidA_IV(α) variant (Table S1) was tested by crossing 25 virgin males with 50 virgin females. This way, 20 different crosses were performed between these seven different lines (Table S2).

Cloning and Sanger sequencing of cidA and cidB variants

The cidA and cidB genes of seven isofemale lines (Luna1, Luna3, Luna8, Luna27, Michele26, Michele1 and Mezzo9) were cloned and Sanger sequenced, as described by Bonneau, Atyame, et al. (2018), starting from the same DNA samples used to determine Wolbachia phylogenetic group. For each gene of each isofemale line, 24 clones were sequenced on average (the detail of numbers of clones sequenced per isofemale line and gene are presented in the Data S2). Moreover, we confirmed the presence of the variants detected in the clones by Sanger sequencing the cidA and cidB fragment amplified from each isofemale line before cloning. This allowed us to verify the polymorphism found in the different cidA and cidB clones. However, even with this double‐checking system, we cannot exclude that some variants might not have been reported. Michele26 was chosen for cloning and sequencing because it was the only incongruent isofemale line found and the six other lines amenable to sustainable rearing under our utilized laboratory conditions. The Muscle alignment tool (Edgar, 2004) implemented in seaview 6.4.1 software (Gouy, Guindon, & Gascuel, 2010) was used to align variant sequences.

The cidB_IV(a/2) variant previously found associated with the incompatible mod phenotype in natural populations from North Africa in the study by Bonneau, Atyame, et al. (2018) was undoubtedly identified only by the cloning and Sanger sequencing because the PCR‐RFLP test (see below) was designed to only discriminate the downstream polymorphic region of cidB variant as it was the one found associated with mod phenotype variation in Bonneau, Atyame, et al. (2018). We therefore named a variant as cidB_IV(a/2) only in situations in which cloning and Sanger sequencing experiments were performed.

Screen of cidA and cidB variants in natural populations from North Italy

Detection of cidA_IV(δ) and cidB_IV(2) variants

We investigated the presence of these variants in the 67 isofemale lines, using the same DNA samples used to determine Wolbachia phylogenetic group. We used the PCR‐RFLP tests described by Bonneau, Atyame, et al. (2018). The cidA_IV and cidB_IV variants described by Bonneau, Atyame, et al. (2018) have two polymorphic regions: an upstream and a downstream region. For instance, cidA_IV(δ/1) and cidA_IV(δ/2) have the same upstream sequence (δ) but two different downstream sequences (1/2), whereas cidB_IV(a/2) and cidB_IV(b/2) have two different upstream sequences (a/b) but the same downstream sequence (2) (Figures S1 and S2). Only the upstream polymorphic region of cidA_IV variants was previously found associated with the CI mod phenotype (Bonneau, Atyame, et al., 2018). The cidA_IV PCR‐RFLP test was, therefore, designed to distinguish between the various upstream polymorphic sequences – cidA_IV(α), cidA_IV(β) and cidA_IV(δ) – regardless of the downstream sequences present. The detection with this test of cidA_IV(δ) in an isofemale line accounts for the presence of cidA_IV(δ/1) and/or cidA_IV(δ/2) (Figure S1). Only the downstream polymorphic region of cidB_IV variants was previously found associated with the CI phenotype (Bonneau, Atyame, et al., 2018). A PCR‐RFLP test was, therefore, designed to distinguish between the cidB_IV(1), cidB_IV(2), and cidB_IV(3) sequences, regardless of the upstream sequences present. Thus, the detection, in isofemale lines, of cidB_IV(2), accounts for the presence of the cidB_IV(a/2) and/or cidB_IV(b/2) variants (Figure S2).

These tests only allowed us to detect variants previously described with the cloning and Sanger sequencing experiment and any variants that were not uncovered with this method would be missed.

In addition to our PCR‐RFLP test, a specific presence/absence PCR test was also designed to detect the presence of cidB_IV(2) (which accounts for the cidB_IV(a/2) and/or cidB_IV(b/2) variants), to confirm the PCR‐RFLP test results in isofemale lines. A 107 bp fragment was amplified with the primers CidB_QPCR_spe_2_dir3 (5′‐GGG‐AAT‐AGT‐GCT‐TTT‐GAT‐AGA‐GAG‐TA) and CidB_QPCR_spe_rev1 (5′‐GTT‐AAA‐CAT‐CTT‐AAA‐CCC‐TCA‐TCA‐CC), under the following PCR conditions: 5 min at 94°C, followed by 35 cycles of 94°C for 30 s, 58°C for 45 s, and 72°C for 30 s, with a final elongation for 5 min at 72°C. To check the specificity, a PCR was performed with CidB_QPCR_spe_2_dir3 and CidB_QPCR_spe_2_rev1 primers on clones carrying either cidB_IV(1), cidB_IV(2) or cidB_IV(3) with the conditions described above and only the clones carrying cidB_IV(2) were amplified demonstrating that our PCR was specific of this variant.

Real‐time quantitative PCR

Quantification of Wolbachia density in male testes

We quantified the density of Wolbachia in the testes of males from the Michele26 and Mezzo9 lines, by quantitative PCR with the LightCycler 480 system (Roche). Specific primer pairs and procedures were described by Berticat, Rousset, Raymond, Berthomieu, and Weill (2002). Each DNA template was obtained from pools of three pairs of testes from six‐day‐old males, with the DNeasy Blood & Tissue Spin‐Column kit (Qiagen; bench protocol: animal tissues). Five independent DNA templates were used for each line (Data S3). We estimated the number of Wolbachia bacteria per mosquito testis, by amplifying two different genes for each sample: the C. pipiens‐specific ace‐2 locus (Weill, Berticat, Raymond, & Chevillon, 2000) and the Wolbachia‐specific single‐copy wsp locus (Berticat et al., 2002). Standard curves were generated with dilutions of a pBluescriptKS vector containing single copies of the ace‐2 and wsp genes. Each DNA template was analyzed in triplicate for the quantification of both wsp and ace‐2. If a triplicate had an error above 0.5 it was removed from the wsp/ace‐2 estimation. As both genes were present as single copies per haploid genome, the ratio of the wsp and ace‐2 signals could be used to estimate the relative number of Wolbachia genomes per Culex genome, thus correcting for mosquito size and DNA quality.

Expression of the cidA and cidB genes

For the Mezzo9 and Michele26 lines, we extracted RNA from 10 six‐day‐old males with Trizol (Life Technologies). The RNA was then treated with DNase with the TURBO DNA‐free Kit (Life Technologies), in accordance with the manufacturer's instructions. The absence of residual DNA was confirmed by performing a PCR specific for cidA and cidB loci with the primers describe in Bonneau, Atyame, et al. (2018). We subjected 2–5 µg of each total RNA sample to reverse transcription with the SuperScript III Reverse Transcriptase Kit and 30 ng of random oligomer primers ([RP]10; Invitrogen, Life Technologies). Four different quantitative PCRs were performed on the resulting cDNA, according to the procedure described by Berticat et al. (2002). The first was specific for the wsp locus, as described by Berticat et al. (2002) and was chosen because (a) it is present in a single copy in the wPipPel reference genome and (b) was the reference gene used for Wolbachia density estimation (Berticat et al., 2002). The second was specific for a 189 bp fragment of the cidA gene conserved in all sequenced wPip strains, and was performed with the primers wPip_0282_QPCR_2_Dir (5′‐AGG‐TCC‐TGT‐ATT‐TGA‐TTT‐CTG‐GA) and wPip_0282_QPCR_2_Rev (5′‐TGA‐ACG‐CGA‐GAA‐AGA‐GCA‐AG). The third was specific for a 135 bp fragment of the cidB gene conserved in all sequenced wPip strains and was performed with the primers wPip_0283_QPCR_1_Dir (5′‐TGA‐GTG‐TTT‐GGA‐GAA‐TGA‐AGG‐A) and wPip_0283_QPCR_1_Rev (5′‐TTC‐CCA‐AAA‐GCA‐AAA‐CCA‐GTT). The fourth was specific for the 107 bp fragment of cidB_IV(2) described above. As wPip strains carry multiple cidB variants we checked that the real‐time quantitative PCR was specific of the cidB_IV(2) variant by performing real‐time quantitative PCR (a) on isofemale lines infected with wPipIV strains lacking cidB_IV(2) but carrying cidB_IV(1) and cidB_IV(3) (Ichkeul 13 and Harash lines) which represent our negative controls that tested for nonspecific amplification and (b) on isofemale lines infected with wPipIV strains carrying cidB_IV(2) (Istanbul and Ichkeul 09) which represented our positive controls. Amplifications were only observed in Istanbul and Ichkeul 09 samples as well as Michele26 and Mezzo9 samples, with melting curve of these samples checked for single product amplification. Each DNA template was analyzed in triplicate for wsp, cidA, cidB and cidB_IV(2). Standard curves were generated for the cidA, cidB, cidB_IV(2) and wsp genes, by diluting the PCR products for these four genes. Expression levels for the cidA, cidB and cidB_IV(2) genes were estimated relative to that of the wsp gene (Data S4).

Statistical analysis

All analyses were performed with r version 3.4.4 software (R Core Team, 2018). Comparisons between the real‐time quantitative variables of the Michele26 and Mezzo9 isofemale lines were performed with the nonparametric Wilcoxon rank‐sum test (Bauer, 1972). The test was chosen because we were comparing two sets of independent data not all normally distributed.

RESULTS

Only cidB_IV(2) variants are associated with the incompatible mod phenotype in North Italy

Both compatible and incompatible mod phenotypes are present in North Italy

A total of 67 isofemale lines were established from larvae collected at four sites in the province of Trento in the North‐East of Italy (San Michele all'Adige, Roveré Della Luna, Mezzocorona and Zambana: Data S1), because it was already known that wPip strains from the wPipIV group occurred in this area (Dumas et al., 2013). The PCR‐RFLP test as described in Altinli et al. (2018) confirmed that all isofemale lines were infected with wPipIV strains. Crosses between males from these 67 isofemale lines and females from the wPipI Tunis laboratory line (reference line used for the screening) led to sort the lines according to the two mod phenotypes previously described (Atyame et al., 2015; Bonneau, Atyame, et al., 2018): the males from 62 isofemale lines from North Italy were found incompatible with Tunis females (qualified as incompatible isofemale lines) and five lines exhibited males compatible with Tunis females (qualified as compatible isofemale lines: Data S1). In summary, 92.5% of the isofemale lines in North Italy exhibited the incompatible mod phenotype.

cidA and cidB variant repertoires

For investigation of the diversity of cidA_IV and cidB_IV genes in North Italian wPipIV‐infected C. pipiens populations and identification of the cidA_IV and cidB_IV variants putatively associated with compatible or incompatible mod phenotypes, we first cloned and Sanger sequenced PCR amplification of cidA and cidB genes from two compatible (Luna 8 and Luna 27) and four incompatible (Luna 1, Luna 3, Michele 1 and Mezzo 9) isofemale lines (Table 1). For each of the six wPip strains studied, we detected several combinations of cidA and cidB variants further referred as repertoires of cidA and cidB variants. To name the different cidA and cidB variants, we used the following nomenclature: the first number corresponds to the group of wPip. Here, as all strains belonged to wPipIV group, all variants were named cidA_IV or cidB_IV. Furthermore in all cidA_IV and cidB_IV variants, two polymorphic regions were detected: the upstream region which is identified with a letter (Greek for cidA variants and Latin for cidB variants) and the downstream region which is identified with a number (Figure S1 [Bonneau, Atyame, et al., 2018]). Consequently, the variants cidA_IV(α/1) and cidA_IV(α/2) share the same sequence for the upstream region but carry a different sequence for the downstream region. The same reasoning was applied to cidB variants.

Table 1

cidA and cidB variant repertoires for seven wPipIV strains from North Italy

Line name(a)	mod phenotype	cidA_IV
		α1	α2	β1	β2	γ1	γ2	δ1	δ2	ε1†	ε2†
Luna 8	Compatible	P	P	A	A	A	A	A	A	P	P
Luna 27	Compatible	P	P	A	A	A	A	A	A	P	P
Luna 1	Incompatible	P	A	A	A	A	A	A	A	P	P
Luna 3	Incompatible	P	A	A	A	A	A	A	A	P	P
Michele 1	Incompatible	P	P	A	A	A	A	A	A	P	P
Mezzo 9	Incompatible	P	P	A	A	A	A	A	A	P	P
Michele 26	Compatible	P	P	A	A	A	A	A	A	P	P

All cidA_IV (a) and cidB_IV (b) variants from Bonneau, Atyame, et al. (2018) and the present study (indicated with †) are compiled. Variants were either present in a wPip strain (P in green) or absent (A in black). Variants matching the mod phenotype are highlighted in bold letters and have a larger font. Michele 26 is in a darker color because it is the incongruent line.

The presence of cidA_IV(α) and the absence of cidA_IV(δ) was observed in the six isofemale lines analyzed (Tables 1a and 2). A new variant called cidA_IV(ε), different from the cidA_IV(β) previously described in North Africa due to the replacement of a valine with an isoleucine residue in position 143, was detected in the six isofemale lines studied (Table 1a; Figure S1). Among all the cidA variants detected in the North Italy wPipIV strains, no specific cidA variant or combination of cidA variants were found associated with either compatible or incompatible mod phenotypes.

Table 2

Distribution of cidA_IV(α), cidA_IV(δ) and cidB_IV(2) in the 247 isofemale lines infected with wPipIV strains from North Africa, Turkey, China and North Italy

cidA cidB variants	mod phenotype
	« Incompatible» n = 79	« Compatible» n = 168
cidA_IV(α)	China: 1/1 North Africa: 15/15 North Italy: 62/62 Turkey: 1/1 Total: 100%	‐ North Africa: 163/163 North Italy: 5/5 ‐ Total: 100%
cidA_IV(δ)	China: 1/1 North Africa: 15/15 North Italy: 0/62 Turkey: 1/1 Total: 21.5%	‐ North Africa: 16/163 North Italy: 0/5 ‐ Total: 9.5%
cidB_IV(2)	China: 1/1 North Africa: 15/15 North Italy: 62/62 Turkey: 1/1 Total: 100%	‐ North Africa: 8/163 North Italy: 1/5 ‐ Total: 5.4%

The prevalence of each variant of interest is given for compatible and incompatible isofemale lines in each area, together with the total percentage of compatible or incompatible isofemale lines carrying the variant concerned, regardless of geographic origin. ‐, indicates the absence of isofemale lines with the compatible mod phenotype in a given location. The isofemale lines from North Africa, Turkey and China were previously analyzed in Bonneau, Atyame, et al. (2018).

‐

‐

‐

‐

‐

‐

The sequencing of cidB gene repertoires of the six wPip strains revealed the presence of cidB_IV(a/2) and cidB_IV(b/2) in the four incompatible isofemale lines and their absence in the two compatible lines (Table 1b). No other cidB_IV variants were found putatively associated with difference in mod phenotypes (Tables 1b and 2). CidB_IV(a/3) and cidB_IV(b/3) were detected in both some compatible and incompatible strains (Table 1b).

For extension of cidA and cidB repertoire analyses at the population scale, we performed a PCR‐RFLP screen of cidA_IV and cidB_IV variants on the 67 isofemale lines originated from the four sites (Data S1). In North Africa natural populations, the upstream region of cidA and the downstream region of cidB were found associated with mod phenotypes variations. Consequently, the PCR‐RFLP tests were designed to differentiate between the different upstream sequences of cidA (α, β, γ, δ) and the different downstream sequences of cidB (1, 2, 3), respectively. Thus, the detection, with this test, of cidA_IV(δ) in an isofemale line accounts for the presence of cidA_IV(δ/1) and/or cidA_IV(δ/2) variants and the detection of cidB_IV(2), accounts for the presence of the cidB_IV(a/2) and/or cidB_IV(b/2) variants. The cidA_IV(α) variant was detected in all the 67 lines. By contrast, the cidA_IV(δ) variant which had been previously reported to be associated with the incompatible mod phenotype (Bonneau, Atyame, et al., 2018) was not found in any of the 67 lines, including the 62 incompatible isofemale lines. The global distribution of the cidA variants of interest in all the natural populations studied in the present study and in those studied in Bonneau, Atyame, et al. (2018) revealed the presence of cidA_IV(α) in all the 247 wPipIV‐infected isofemale lines studied, regardless of mod phenotype and geographic origins (Table 2; Figure 1). Overall, cidA_IV(δ) was detected in 21.5% (17/79) of the incompatible isofemale lines (Table 2; Figure 1). All together these data show a lack of correlation between cidA variants and mod phenotype variation.

All the 62 incompatible Italian isofemale lines carried the cidB_IV(2) variants, as confirmed by the two independent PCR‐based methods (Data S1). However, this variant was also detected in Michele26, one of the five compatible isofemale lines. In compiling the data from China, Turkey, North Africa and North Italy, cidB_IV(2) was detected in 100% (79/79) of the incompatible isofemale lines while only in 5.4% (9/168) of the compatible isofemale lines, regardless of geographic origin (Table 2).

All lines that carry cidA_IV(α) are reciprocally compatible

Because the cidA_IV(α) was detected in all studied wPipIV strains and all tested mosquito lines infected with wPipIV strains were always found mutually compatible (Atyame et al., 2014), we previously hypothesised that this variant might be involved in this reciprocal compatibility and thus in the resc function (Bonneau, Atyame, et al., 2018). To investigate further that hypothesis, crosses between all the different mosquito wPipIV infected lines present in the laboratory were performed. We were able to perform a total of 20 new crosses between wPipIV‐infected isofemale lines originated from Turkey, North Africa and Italy which carry the cidA_IV(α) variant. These 20 crosses together with the 38 already performed in Atyame et al. (2014) confirmed the compatibility between wPipIV strains (Table S2).

How can Michele26 be compatible while it carries the cidB_IV(a/2) variant associated with incompatibility?

The Michele26 isofemale line harbours a unique cidB variant repertoire

The cidA and cidB variant repertoires of the wPipIV strain infecting the Michele26 isofemale line were cloned, Sanger sequenced, and compared with the repertoires obtained for the four incompatible and two compatible isofemale lines from North Italy (Table 1a). No specific cidA_IV variant repertoire was identified for this isofemale line. Indeed, the cidA_IV variant repertoire was identical to that of some compatible and incompatible isofemale lines (Table 1a). In contrast the cidB_IV variant repertoire of Michele26 was unique by including cidB_IV(a/2), but lacking cidB_IV(b/2), which was present in all the incompatible lines. Furthermore, two variants unique to Michele26, so called cidB_IV(c/1) and cidB_IV(c/3), were identified (Table 1b; Figure S2). In summary, the Michele26 isofemale line had a similar cidA variant repertoire but a unique cidB variant repertoire different from the other compatible and incompatible wPipIV strains sampled in the same area.

Michele26 males are able to induce CI

The incapacity of males Michele26 to induce CI when crossed with Tunis females might be due to the incapacity of Michele26 males to induce CI at all. Consequently, we checked the capacity of Michele26 males to induce CI by crossing them with females artificially cured from their Wolbachia by tetracycline treatment. A total of 21 eggs‐rafts were collected and none of them hatched demonstrating that Michele26 males were able to induce CI.

cidB_IV(2) expression is lower in compatible Michele26 males than in incompatible Mezzo9 males

We then investigated possible differences in the expression of the cidA and cidB genes between compatible Michele26 males and incompatible Mezzo9 males. Indeed, the compatible phenotype of Michele26 could result from an absence of expression of the cidB_IV(2) variant associated with incompatibility. The overall levels of cidA and cidB expression between the two lines were not significantly different (cidA Wilcoxon, W = 62, p = .182 and cidB Wilcoxon, W = 60, p = .243, Table 3; Data S4). The expression of cidB_IV(2) was studied by real‐time quantitative PCR of a sequence fragment accounting for both the cidB_IV(a/2) and cidB_IV(b/2) variants (a/2 present in Mezzo9 and Michele26, and b/2 present only in Mezzo9, see Table 1b). Expression of the cidB_IV(2) fragment in males from the Michele26 isofemale line was significantly lower than in males from the Mezzo9 isofemale line (0.04 as opposed to 0.06) (Wilcoxon, W = 12, p = .003, Table 3; Data S4).

Table 3

Total cidA and cidB variant expression levels, specific cidB_IV(2) expression levels and Wolbachia density in the testes of compatible Michele26 and incompatible Mezzo9 males

Line	Ratio
	Wolbachia testes density (wsp/ace−2)	cidA expression (cidA/wsp)	cidB expression (cidB/wsp)	cidB_IV(2) expression (cidB_IV(2)/wsp)
Mezzo9	18.10 ± 4.72 (a)	1.67 ± 0.56 (a)	0.72 ± 0.23 (a)	0.06 ± 0.01 (a)
Michele26	10.25 ± 3.42 (b)	1.33 ± 0.31 (a)	0.58 ± 0.18 (a)	0.04 ± 0.02 (b)

For each variable, letters indicate groups of statistical similarity. Mean values are expressed ± standard deviation. Wolbachia density and expression level were estimated using real‐time quantitative PCR as ratio between target genes (wsp for Wolbachia density, cidA, cidB and cidB_IV(2) for expression level) and reference genes (ace‐2 for Wolbachia density and wsp for expression level).

Less Wolbachia in the testes of males from the compatible Michele26 isofemale line

As CidB proteins are predicted to be introduced in the sperm during spermatogenesis, we determined Wolbachia density in the gonads of both Michele26 and Mezzo9 males. The testes of males from the Michele26 isofemale line contained significantly less Wolbachia than those of males from the Mezzo9 isofemale line: 10.25 ± 3.42 and 18.10 ± 4.72 Wolbachia per host cell, respectively (Wilcoxon, W = 1, p = .008, Table 3; Data S3).

DISCUSSION

In the current state of knowledge on CI and its diversity in C. pipiens, new investigations on the putative role(s) of cidA Pip were necessary. Indeed, the fact that the coexpression of both cidA Pip and cidB Pip in D. melanogaster males was required to induce CI could support the implication of both cidA Pip and cidB Pip in the mod function (Figure 1; LePage et al., 2017; Shropshire & Bordenstein, 2019; Shropshire et al., 2018). However, the same requirement for the production of live transgenic D. melanogaster and S. cerevisiae could suggest that CidA Pip protein may simply serve as an antidote (i.e., resc) to CidB Pip protein without being directly involved in the mod function (Figure 1; Beckmann et al., 2019b, 2017). As we could not conduct a functional transgenic study of the role of cidA in C. pipiens, due to technical restrictions and the amplifications of this gene in the Wolbachia harboured by this species, we investigated here the putative link between CI mod phenotype diversity and variation in its cidA repertoire by sampling new natural populations in North Italy infected with wPipIV Wolbachia, a wPip phylogenetic group for which simple mod phenotype variations were already screened in North Africa (Atyame et al., 2015; Bonneau, Atyame, et al., 2018). The screening of the 67 isofemale lines obtained from our sampling in Italy revealed the coexistence of the two mod phenotypes in these Italian populations. However, unlike North African populations in which 8.4% of the isofemale lines were found incompatible when males from these lines were crossed with wPipI‐infected females from the Tunis line (Atyame et al., 2015; Bonneau, Atyame, et al., 2018), 92.5% of the isofemale lines from North Italian populations were found incompatible.

In natural populations from North Africa, both cidA_IV(δ) and cidB_IV(a/2) variants were associated with the incompatible mod phenotype (Table 2; Figure 1), a pattern that suggest that both cidA and cidB, were putatively involved in the mod phenotypic variations (Bonneau, Atyame, et al., 2018). In North Italy, cidB_IV(2), which accounts for both the cidB_IV(a/2) and cidB_IV(b/2) variants, was also found associated with the incompatible mod phenotype (Tables 1b and 2; Figure 1; Data S1). All our data demonstrate that only cidB_IV(2) variants were systematically found in incompatible repertoires, strengthening the link between cidB variations and mod phenotype variations in C. pipiens. By contrast to the natural populations from North Africa, cidA_IV(δ) was not detected in any of the wPipIV strains hosting by the 67 isofemale lines in North Italy, demonstrating that this variant was not essential for the incompatible mod phenotype (Tables 1a and 2; Figure 1). We were unable to identify any other cidA variant or combination of variants associated with either incompatible or compatible mod phenotypes (Table 2; Figure 1). Furthermore, we found exactly the same cidA repertoire associated with either incompatible or compatible mod phenotypes, suggesting that cidA plays no role in mod phenotype diversity in C. pipiens (Tables 1a and 2). We can thus hypothesize that the association between cidA_IV(δ)/cidB_IV(2) and the incompatible mod phenotype in natural populations from North Africa resulted from codiversification of these two variants. As previously suggested, cidA and cidB may encode a toxin‐antidote (TA) system in which CidA acts as the antidote of CidB (Beckmann et al., 2019b, 2019a, 2017; Shropshire et al., 2019). Such TA system may have driven the association of cidA_IV(δ) and cidB_IV(a/2) variants in North Africa if these variants interact particularly well together independently of any involvement of cidA in the mod phenotype diversity. However, in the absence of cidA_IV(δ) in North Italian isofemale lines, other cidA_IV variants may also interact with cidB_IV(2).

In all the 247 wPipIV infected C. pipiens lines yet investigated (Table 2), including the 67 North Italian ones, the cidA_IV(α) variant was detected (Figure 1). A total of 58 intra‐wPipIV group compatible crosses between Turkish, North African and Italian lines, including the 20 crosses from the present study, show the self‐compatibility between wPipIV‐infected isofemale lines previously established (Atyame et al., 2014). As already suggested in Bonneau, Atyame, et al. (2018), this observation supports a role for cidA in the resc function in C. pipiens, as the presence of a ubiquitous cidA variant is expected to explain the compatibility of mosquitoes infected with wPipIV strains. This conclusion is further supported by the recent findings of Shropshire et al. (2018), revealing the involvement of cidA Mel in the resc function in transgenic D. melanogaster females (Table 2).

Our results show that cidB Pip variant repertoire is associated with the diversity of mod phenotypes observed in C. pipiens. Together with functional transgenic data (Beckmann et al., 2017), they clearly demonstrate the involvement of cidB in both mod function and mod phenotype diversity in C. pipiens (Figure 1). By contrast, we show here that cidA Pip is not involved in mod phenotype diversity, as lines with different mod phenotypes had the same cidA Pip repertoire. In C. pipiens we have thus far no proof of a two‐by‐one system as proposed for D. melanogaster (Shropshire & Bordenstein, 2019; Shropshire et al., 2018). The cidA and cidB genes may not, therefore, behave in the same way in the Wolbachia bacteria infecting C. pipiens and D. melanogaster. Further investigations of this putative divergence in the molecular mechanisms of CI induced by cidA Mel and cidA Pip in these two species are required to shed light on this putative difference.

In North Africa, ~5% of isofemale lines were incongruent, i.e., exhibiting a compatible mod phenotype while carrying the cidB_IV(a/2) variant associated with incompatible mod phenotype. This phenomenon could not be further studied in our previous work in North Africa as the lines were not alive anymore when the variants were screened. In Italy, we successfully sampled and maintained one incongruent isofemale line (Michele26). We investigated this discordant isofemale line further to search for possible causes of this apparent dissociation between genotype and phenotype. First of all, the incapacity of Michele26 males to induce CI when crossed with Tunis females could have resulted from the incapacity of males Michele26 to induce CI at all. However, we confirmed that Michele26 was able to induce CI by crossing males with females artificially cured of Wolbachia from the SlabTC laboratory line. More importantly, we showed that the wPipIV strain harboured by Michele26 mosquitoes presented specific genetic features distinguishing it from both compatible and incompatible lines found in Italy. We found a specific cidB variant repertoire, including two variants (cidB_IV[c/1] and cidB_IV[c/3]) not detected in any other wPipIV strains cloned and Sanger sequenced (Table 1b; Figure S2). Such presence of additional cidB_IV variants in Michele26 might result in the deregulation of titration or binding with other Wolbachia or host targets, preventing incompatibility with wPipI females. The wPipIV strain harboured by Michele26 also lacked the cidB_IV(b/2) variant reported in all incompatible isofemale lines from North Italy. It was tempting to speculate that both cidB_IV(a/2) and cidB_IV(b/2) were required to induce incompatible phenotype. However, this hypothesis was ruled out by the lack of detection of cidB_IV(b/2) in incompatible isofemale lines from North African populations. Michele26 males do express cidB_IV(2), so that compatibility cannot be caused by a lack of expression of this variant. However, the levels of cidB_IV(2) expression in compatible Michele26 males were lower compared to those in incompatible Mezzo9 males. This is certainly because Michele26 carries only the cidB_IV(a/2) variant, whereas Mezzo9 carries both cidB_IV(a/2) and cidB_IV(b/2) variants. As CidB proteins are probably released into the sperm in the testes during spermatogenesis, we measured the density of Wolbachia in the male gonads. We found that the density of Wolbachia in Michele26 was significantly lower than in Mezzo9. The lower cidB expression level in addition with a lower density of Wolbachia in males Michele26 could result in an insufficient amounts of CidB_IV(2) proteins being produced to induce incompatibility with wPipI females. A similar dosage‐driven hypothesis had already been proposed in a quantitative model where CI phenotype diversity could rely on different mod and resc genes as well as the amount of these genes products (Nor et al., 2013).

In conclusion, our findings support that cidB Pip variant repertoire is associated with the diversity of mod phenotypes observed in C. pipiens. Together with functional transgenic data (Beckmann et al., 2017), they clearly suggest the involvement of cidB in both mod function and mod phenotype diversity in C. pipiens but further suggest that variation in wPip density and/or cidB expression may matter. By contrast, we have no indication that cidA Pip could be involved in mod phenotype diversity, as lines with different mod phenotypes exhibited exactly the same cidA Pip repertoire. Overall, a toxin‐antidote model where cidB is a toxin and cidA its antidote fits well our current knowledge of C. pipiens‐Wolbachia interactions.

AUTHOR CONTRIBUTIONS

M.B., M.S., and M.W. conceptualized and designed the study; M.B., B.C., D.A., M.S., and M.W. field sampled the mosquitoes; M.B., A.L., M. P‐S, R.C., S.U. performed the experiments; M.B., M.S., and M.W. analyzed and interpreted the data; M.B., B.C., R.C., M.S., and M.W. wrote the manuscript.

Click here for additional data file.

Click here for additional data file.

Click here for additional data file.

Click here for additional data file.

Click here for additional data file.

Click here for additional data file.

Click here for additional data file.

Click here for additional data file.

Click here for additional data file.

Click here for additional data file.