Genetic Diversity of the Invasive Gall Wasp Leptocybe invasa (Hymenoptera: Eulophidae) and of its Rickettsia Endosymbiont, and Associated Sex-Ratio Differences

The blue-gum chalcid Leptocybe invasa Fisher & LaSalle (Hymenoptera: Eulophidae) is a gall wasp pest of Eucalyptus species, likely native to Australia. Over the past 15 years it has invaded 39 countries on all continents where eucalypts are grown. The worldwide invasion of the blue gum chalcid was attributed to a single thelytokous morphospecies formally described in 2004. Subsequently, however, males have been recorded in several countries and the sex ratio of field populations has been found to be highly variable in different areas. In order to find an explanation for such sex ratio differences, populations of L. invasa from a broad geographical area were screened for the symbionts currently known as reproductive manipulators, and both wasps and symbionts were genetically characterized using multiple genes. Molecular analyses suggested that L. invasa is in fact a complex of two cryptic species involved in the rapid and efficient spread of the wasp, the first recovered from the Mediterranean region and South America, the latter from China. All screened specimens were infected by endosymbiotic bacteria belonging to the genus Rickettsia. Two closely related Rickettsia strains were found, each infecting one of the two putative cryptic species of L. invasa and associated with different average sex ratios. Rickettsia were found to be localized in the female reproductive tissues and transovarially transmitted, suggesting a possible role of Rickettsia as the causal agent of thelytokous parthenogenesis in L. invasa. Implications for the variation of sex ratio and for the management of L. invasa are discussed.


Introduction
Leptocybe invasa Fisher & La Salle (Hymenoptera: Eulophidae: Tetrastichinae) commonly known as blue-gum chalcid is a gall wasp of many Eucalyptus species. It was unknown until the

Ethics Statement
The sampling of living material involved in our experiments included wasps, i.e. Leptocybe invasa, associated with galls on Eucalyptus sp. All sampling locations are not privately owned or protected (coordinates in Table 1). Besides neither the host plant nor the wasp species are endangered or otherwise protected, and therefore no specific permits were required for these locations/activities.

Insect sampling
Stands of Eucalyptus camaldulensis Dehnhardt or E. globulus Labillardière infested with L. invasa were sampled from Portici, S. Maria al Bagno and Siracusa (Italy), La Plata (Argentina) and Hubian (China). Additional material was received from Antakya (Turkey) and Mouaden (Tunisia) ( Table 1). Leaves harbouring galls of different ages were removed from the trees, placed in plastic boxes and stored at room temperature until the emergence of adults. Soon after their emergence, wasps were killed in 95% ethanol and morphologically identified using the original species description of females [3] and males [21].

Molecular characterization and phylogenetic analyses of L. invasa
DNA was extracted from whole single individuals (males and females) by using both destructive and non-destructive Chelex and proteinase K based methods modified as in [36]. In this latter case, wasps were then treated and card mounted as in [37]. Three genes were sequenced: the mitochondrial cytochrome c oxydase subunit I (COI), and two ribosomal genes, the expansion segment D2 of the 28S ribosomal subunit (28S-D2) and the Internal Transcribed Spacer 2 (ITS2). COI is a very effective and universally used marker for species delimitation [38]. ITS2 has been successfully used at the species level in many Chalcidoidea taxa, e.g. Pteromalidae [39], Trichogrammatidae [40] and Eulophidae [37]. 28S-D2 has also proved diagnostic at the species level in Pteromalidae [39], Aphelinidae [41] and Eulophidae [36]. Two regions of COI were amplified, the first using the forward primer LCO1490 paired with either HCO2198 [42] or Lep-r1 [38] reverse primers, and the second pairing the forward primer C1-J-2183 [43] or Apf [44] with the reverse primer TL2-N-3014 [43]. The ribosomal gene 28S-D2 was amplified with primers D2F and D2R [39]. PCR reactions and cycling conditions for COI and 28S-D2 were set as described in [36]. For the amplification of ITS2, primers ITS2F [39] and ITS2Tri-chRev [40] were used in PCR reactions as in [45]. PCR products were checked on a 1.2% agarose gel stained with ethidium bromide and directly sequenced. Chromatograms were assembled, edited by eye and aligned with Bioedit 7.2.5 [46] and the sequences deposited in GenBank, with the accession numbers reported in Table 1.
All sequences with chromatograms showing ambiguous peaks were cloned. Amplicons were ethanol precipitated, ligated into the pGEM-T Easy plasmid vector (Promega), and cloned into Escherichia coli TOP10 competent cells (Invitrogen) according to the manufacturer's instructions. Transformants were PCR-screened with universal M13 vector primers, and inserts of the expected size were sequenced. COI sequences were virtually translated into the corresponding amino acid chain to detect frame-shift mutations and stop codons, which may indicate the presence of pseudogenes [47,48].
Phylogenies for L. invasa were reconstructed using Bayesian inference (BI), utilizing MrBayes 3.1.2 [49], following the methodology reported in [37] on concatenated (ITS2-28S-D2-COI) alignment by using the GTR+G+I nucleotide model as selected by jModeltest [50]. Sharing some morphological features with L. invasa, two species of closely related genera,    Baryscapus silvestrii Viggiani et Bernardo (Hymenoptera: Eulophidae) and Aprostocetus monacoi Viggiani (Hymenoptera: Eulophidae) were used as outgroups to root the tree. Sequence distances within and between populations of COI dataset were calculated as uncorrected pairwise distances (p-distances) using MEGA 4 [51]. We also looked for diagnostic single nucleotide polymorphisms (SNPs), i.e. character state in a given nucleotide position shared by all individuals from one group and different from all individuals in any other group, in COI sequences, as well as the occurrence of diagnostic non-synonymous amino acid changes. For species delimitation, we used the Poisson tree processes model (PTP) [52], a recently developed method that been successfully applied [53]. All sequences of Tetrastichinae (5'-COI "barcoding" region) available in GenBank were included, aligned in BioEdit and run in MrBayes (GTR+G+I substitution model) for 1 million generations (rate sample = 1000, burn in value = 25%). BI tree obtained was run on the bPTP web server as in [53], for 500,000 MCMC generations. The relationships between L. invasa specimens were also investigated using Statistical Parsimony in TCS 1.21 [54] on the COI dataset, as in [55].

Molecular characterization of Rickettsia
Multiple genes were sequenced to characterize Rickettsia at the strain level. Approximately 1000 bp of the 16S rRNA gene were amplified with Rb-F and Rb-R primers as previously described;~600 bp of the citrate synthase (gltA) gene were obtained by a nested PCR strategy with primers CS1d and CS1273r [64] followed by CS78d paired with CS715r [64]; the subunit α of ATP synthase F1 (atpA) gene (~750 bp) was amplified with ATP syn f1 α fw and ATP syn f1 α rev primers [65]. For both genes, amplifications were carried out using a "step-up" program: 5 min of initial denaturation at 94°C, first 10 cycles step at 94°C for 40 sec, 37°C for 2 min and 72°C for 90 sec and the second 35 cycles step with the same denaturation and elongation parameters and annealing temperature set at 48°C for 1 min. The amplification was completed by holding for 5 min at 72°C. Finally, the intergenic spacer rpmE-tRNAf Met , known to be useful for characterizing strains of Rickettsia [66], was sequenced using rpmEF and rpmER primers [67]. All sequences were assembled, edited and aligned with Bioedit 7.2.5 [46], and deposited in GenBank, with the accession numbers reported in Table 1.

Rickettsia phylogeny and host-endosymbiont relationship
Sequences of the 16S rRNA, gltA and atpA genes were aligned with homologous sequences available in GenBank. The related species Orientia tsutsugamushi (Hayashi) was included as outgroup in the 16S rRNA and atpA genes analyses [68], whereas a Wolbachia strain was used as outgroup in the gltA gene analysis. For Rickettsia 16s rRNA, gltA and atpA genes, phylogenetic reconstructions using Bayesian inference was performed using MrBayes 3.1.2 [49], following the methodology reported in [37]. Intergenic spacer region rpmE-tRNAf Met sequences of L. invasa were compared to each other and checked against GenBank for highest similar matches.
To test for congruence between the host and endosymbiont trees, trees of L. invasa populations and of their own Rickettsia symbionts were compared. For this analysis, ML trees were built including P. soemius and its own Rickettsia sequences as outgroups, as the only eulophid taxon currently available in GenBank for which all five insect and symbiont genes here used are available.

Fluorescence microscopy
Localization of Rickettsia within the host's reproductive tissues was carried out, for the Italian population, using fluorescence in situ hybridization (FISH) with the Rickettsia-specific probe RickPn-Cy3 [33] and the universal probe EUB338-Cy5 [69]. Ovaries and mature eggs were extracted from adult females in a drop of PBS buffer under a stereomicroscope. Ovaries and eggs were subjected to the whole-mount FISH method reported by [33] except for the duration of egg dechorionation (10 min in 50% commercial bleach in PBS solution). Stained samples were observed both under a ZEISS Axiophot 2 epifluorescence microscope and under the Leica confocal TCS SP5 microscope. Images obtained on multiple focal planes were stacked with software Leica application Suite, Advances fluorescence 2.4.1. The specificity of the observed signals was verified with the following control experiments: no-probe control and RNase-digested control. Nuclei of the host cells were counterstained with DAPI (0.4 μg/ml) in mounting medium.

Morphological characterization of Leptocybe invasa
All female and male specimens from different countries were identified as L. invasa, and did not show any appreciable difference in their morphology. Voucher material (listed in Table 1) is deposited at the Institute for Sustainable Plant Protection (IPSP, Portici, Italy).

Molecular characterization and phylogenetic analyses of L. invasa
The ribosomal gene 28S-D2 sequences were obtained for 27 specimens (Table 1). Sequences ranged from 488 to 578 nt, and, after trimming, the final alignment consisted of 482 bp. The ribosomal gene 28S-D2 sequences of Italian, Tunisian, Argentinean and Turkish samples were identical to each other and differed by a single nucleotide from Chinese sequences (Table 2). Similarly, ITS2 sequences from Italian, Tunisian, Argentinean and Turkish samples were identical, whereas the Chinese sequences had a single missing nucleotide (427 bp versus 428 bp) ( Table 2).
Two regions of COI (1455 nucleotides in total) were sequenced from each of 19 specimens ( Table 1). The same COI haplotype was recovered from the Italian, Argentinean and Tunisian populations, whereas one and two unique haplotypes were recovered from the Chinese and Turkish populations, respectively ( Table 2). Clones of COI amplicons, that showed ambiguous peaks in the chromatograms of sequences from the Turkish population, mostly resulted from pseudogenes, except for two clones that were then used in subsequent analyses. BI phylogenetic analyses of the concatenated ITS2, 28S-D2 and COI (S1 Fig) showed two well-supported clades. The first clade includes L. invasa from China, and the second clade includes Italian, Argentinean and Tunisian populations (hereafter called lineage A) and the Turkish population.
The bPTP analysis indicates a moderate support for the existence of two cryptic species in the blue gum chalcid complex, delimiting L. invasa from China as a putative species (hereafter called Chinese lineage) and the monophyletic clade including lineage A and Turkish specimens as a second putative species (hereafter called Western lineage) (Fig 1). Based on this analysis, there is 65% probability that the four COI haplotypes do not represent a single species, 49% probability that the 3 haplotypes "lineage A", "TK1" and "TK2" represent a single species (which is consistent with the slightly high uncorrected intra-specific distances in the Turkish population, see below) and 62% probability that the Chinese population is a separate species. Statistical parsimony on the COI dataset yielded two separate networks, corresponding to the Western and Chinese lineages. The connection limit necessary to obtain a single network was 48 steps.
The distance between Chinese and Turkish samples was 3.7%, that between Chinese lineage and lineage A was 3.1% (Table A in S1 File), and that between Chinese and Western lineages was 3.6% (Table B in S1 File). The two Turkish haplotypes resulted in an intraspecific distance of 1.5% (Table C in S1 File). Analysis of diagnostic sites for COI sequences showed that the Chinese lineage is distinguishable by 43 SNPs (S3 Table). Three diagnostic non-synonymous amino acid changes were found in the Chinese lineage: valine instead of isoleucine, methionine instead of isoleucine and serine instead of threonine.

Detection of bacterial symbionts
DGGE analysis of PCR-amplified 16S rRNA gene showed that the only endosymbiont infecting all tested L. invasa adults (26 females and 4 males) was Rickettsia (data not shown). Congruent results were obtained by PCR screening, which revealed that 100% of individuals were infected only by Rickettsia.

Molecular characterization of Rickettsia
Rickettsia 16S rRNA gene sequences (ranged from 744 to 837 bp) were obtained from 16 specimens (Table 1). No intra-lineage variation was found except within the Turkish population, where two sequences were recovered: one in specimen Li-TK1, shared with lineage A, and another found exclusively in specimen Li-TK2 (different by two nt). 16S rRNA gene of the Chinese lineage differed from all the others by two nt ( Table 2). Sequences of the Rickettsia gltA gene (499 bp), obtained from 18 wasps (Table 1), were invariant, except for in the Chinese lineage, which differed by a single nucleotide. Sequences of the atpA gene (681 bp), obtained from 15 wasps (Table 1), were identical in Rickettsia found in all L. invasa specimens. rpmE-tRNAf-Met sequences (381 bp) were obtained from 11 wasps (Table 1). Rickettsia in individuals of the Western lineage showed the same sequence, which differed from the sequence found in Rickettsia in the individuals of the Chinese lineage by three nucleotides. The Chinese and all other samples sequences had a match of 97 and 98% identity with R. felis rpmE-tRNAf Met sequences, respectively.
Rickettsia phylogeny and host-endosymbiont relationship BI phylogenetic analyses of 16S rRNA, gltA and atpA genes showed that Rickettsia symbionts of L. invasa form a lineage in the Rickettsia transitional group (Fig 2 and S2 Fig).

Fluorescence microscopy
Using FISH analysis, we found that Rickettsia occur in the reproductive tissues of L. invasa females. Inside the ovary, dense clusters of bacteria were observed in the ovarioles ( Fig 3A) and in the germarium (Fig 3C). Simultaneous probing with the Rickettsia specific probe and the universal bacterial probe (S4A and S4B Fig) did not reveal the occurrence of different bacteria out of Rickettsia ( Fig 3B) showing a full overlay of two probes. Within the developed oocytes, bacteria were densely distributed both near the nucleus (Fig 3D and 3E) and in the peduncle (Fig 3F). Negative controls (no-probe and RNase-digested controls) did not display any signal, confirming the specificity of the detected signals.

Discussion
Specimens of L. invasa from different countries here examined all correspond to a single morphospecies. However, the integration of phylogenetic analysis, bPTP species delimitation method, statistical parsimony, COI distances and diagnostic SNPs, genetic differences at the strain level of Rickettsia symbionts, and sex ratio differences supports separation of the Chinese and Western lineages, which can be therefore treated as putative species. Indeed, two genetically distinct lineages of L. invasa were identified on the basis of molecular, phylogenetic and species delimitation analyses. One lineage included populations from the Mediterranean region and South America, the other was present in China. The two lineages showed differences in all sequenced genes. Whereas ITS2 resulted to be more conserved than in other chalcidoid taxa, variation in COI and 28S-D2 is compatible with species level divergence reported previously. Indeed, 28S-D2 is often invariant between closely related species that are instead well differentiated on the base of other genes (usually COI) and biological traits [70,45]. In one case, even Palearctic and Nearctic species of the genus Pnigalio (Eulophidae) shared the same 28S-D2 sequence [45]. The single 28S-D2 polymorphism in L. invasa is consistent with what already observed in Eulophidae, where a high variation in COI can correspond to just a single diagnostic polymorphism in 28S-D2 [36]. As for COI, the divergence (<4%) is similar to that between two species pairs of the genus Encarsia (Aphelinidae) (Gebiola et al. in prep). The lack of COI genetic variation in specimens from the same locality (with the exception of the Turkish population) could be explained either by founder effects (the reduced genetic variation that occurs when a population is established starting from a few or a single specimen) [71], or by endosymbiont infection. Endosymbionts, acting as reproductive manipulators, are considered responsible for the very low mitochondrial genetic diversity in infected populations [32,72], indeed they can induce selective sweeps that indirectly impact the existing polymorphism of mtDNA [73,74].  Among the known reproductive manipulators, Rickettsia was the only one found to infect L. invasa. Phylogenetic reconstructions based on three different genes (16S rRNA, gltA and atpA) placed this symbiont within the Rickettsia transitional group, in a clade that includes the symbionts of the eulophid wasps Pediobius rotundatus (Fonscolombe) and N. formosa [68]. Congruently, the rpmE-tRNAf Met gene sequence showed the highest similarity with R. felis, a species also belonging to the transitional group [68,75]. In particular, Rickettsia gene sequence analysis revealed the occurrence of two closely related bacterial strains, one associated with the putative Western species and one with the putative Chinese species. The 16S rDNA, gltA, atpA and rpmE-tRNAf Met sequences of the symbionts of the Western species samples showed no differences with the exception of 16S rRNA gene sequence of the Turkish specimen Li-TK2, differing by two nt. Instead, the Rickettsia strain associated with the Chinese species differed from the Western species in three out of four genes (16S rRNA, gltA and rpmE-tRNAf Met ).
The strict association between the two Rickettsia strains and the two L. invasa lineages is also evident by the striking congruence between symbiont and host phylogenetic trees (S3 Fig), which suggests strict vertical transmission of the symbiont. Congruent topologies of host and symbiont molecular phylogenies are consistent with coevolution, but could be caused by other factors as well, such as resource tracking, where the symbiont evolves in response to a host trait with a phylogenetic signal [76].
Populations of L. invasa are known to have different sex ratios over their geographic distribution. Males have never been recorded from Italy ( [1], see also S2 Table), Tunisia [34] and Argentina [35], whereas males are rare in Turkey (sex ratio 0-0.5%) [21]. In these countries, it is clear that L. invasa reproduces by thelytokous parthenogenesis. In contrast, a high frequency of males has been recorded in populations from China (males ranging from 18-48% in some populations) [17,26] and in populations from South-East Asia (Taiwan, India, Thailand) [22,23,24,25]. Interestingly, the genetic differentiation within L. invasa shown here mirrors the geographic variation in sex ratio. However, despite the occurrence of males, there is evidence suggesting that some populations from South-East Asia also reproduce by thelytokous parthenogenesis (although biparental populations cannot be excluded). Indeed, populations with only occasional males occur in China as well [77]. Furthermore, virgin females of the Chinese population studied in this work (Wang, unpublished observations) and of a Thai population [24] are able to produce both females and males, with male progeny apparently being not functional [24].
Why do such differences in the sex ratio exist between thelytokous populations of the Western and Chinese putative species? Results discussed above and the evidence that Rickettsia symbionts of L. invasa occur at high density within the ovaries and the eggs and is transmitted vertically to the progeny with very high efficiency support the hypothesis that Rickettsia can induces thelytokous parthenogenesis in L. invasa. Removing Rickettsia symbionts by feeding females antibiotics should restore the production of male progeny [33,78], but, despite numerous attempts, we were not able to rear and thus cure infected L. invasa. Furthermore, the interaction between host and symbiont genotype plays an important role in the phenotypic effect and transmission efficiency of reproductive manipulators [79]. Low transmission efficiency of Rickettsia could make the sex ratio of thelytokous individuals more susceptible to the effect of environmental factors. Efficiency of symbiont transmission through the host germline, but also penetrance of the reproductive phenotype, and infection prevalence in the host population, may be correlated with bacterial density [80]. Bacterial density is in turn regulated by genetic factors of the host and the symbiont itself and is strongly influenced by environmental factors, like temperature, antibiotics, and host age [81,82]. The general variation in bacterial density in response to temperatures indicates that there can be large spatial, temporal, and seasonal differences in endosymbiont densities and functions in natural populations. High temperatures, by reducing the symbiont density within the reproductive tissues of the host, can induce the production of male progeny by infected parthenogenetic females and variable sex ratios in field populations [83,84]. In the case of L. invasa, laboratory experiments with a Chinese population have shown that at constant temperatures the frequency of males in the progeny of thelytokous females increases from 2% at 20-23°C to 7% around 30°C [77]. Therefore, it is possible that the association between Rickettsia strain and its Chinese host is weaker (in terms of bacterial density and/or transmission efficiency) than that occurring with its Western host, and consequently more susceptible to the effect of high and/or low temperatures. Another plausible explanation for the different sex ratios of the Chinese species may be a more recent origin of the symbiotic association with Rickettsia than in the Western species. Thus, the production of frequent male offspring by infected Chinese females may be due to a maladaptive side effect of incomplete coevolution between host and symbiont, as recorded in another host-parasitoid system [85]. Besides, a possible influence of the host plant cannot be excluded, as the Chinese specimens were collected on E. globulus instead of E. camaldulensis. It has been demonstrated in other systems that plants on which phytophagous insects feed may influence the host-symbiont relationship [86,87]. Lastly, unless a direct involvement of Rickettsia in the thelytokous reproduction [33,78] or in the oogenesis [88] of L. invasa is conclusively demonstrated (by obtaining males or no eggs respectively from cured females), we cannot rule out the possibility that thelytoky could be genetically determined.
The likely existence of more than one species behind the morphospecies L. invasa could have important implications in terms of pest management. For example, parasitoids introduced in several countries [19] could have a different degree of specificity towards the two cryptic species, being able to parasitize only one of them or performing sub-optimally on different host species. It is therefore very important to determine whether the host range of the main parasitoid species described as monophagous (Selitrichodes kryceri Kim & La Salle, Selitrichodes neseri Kelly & La Salle, and Quadrastichus mendeli Kim & La Salle) [4,89] includes both cryptic species of L. invasa. Moreover, an influence on parasitoid performance might be caused not just by different hosts but also by the presence of diverse bacteria, as recently showed in other systems where different species of endosymbionts, or even slightly different strains, have an impact on the specificity of host-parasitoid interactions [90][91][92][93].

Concluding Remarks
We showed that the world wide distribution of L. invasa very likely involves at least two species showing distinct sex-ratios, that Rickettsia may be the causal agent of thelytokous reproduction, that two different symbiont strains are associated with the two putative host species, and that the host evolutionary history is recapitulated in the relationship of their microbial symbionts. Based on these results, a better understanding of the interaction between the host and symbiont is critical to explain biological differences. Furthermore, a deeper knowledge and characterization of the different populations of L. invasa from around the world is an essential challenge that should be addressed because of its consequences on pest management. As L. invasa is now widespread and biological control seems to be the best solution of its management, it is important to reassess the efficiency of the parasitoids currently used on both cryptic gall wasps, in order to avoid failures of biological control programs.