Ooctonus vulgatus (Hymenoptera, Mymaridae), a potential biocontrol agent to reduce populations of Philaenus spumarius (Hemiptera, Aphrophoridae) the main vector of Xylella fastidiosa in Europe

Sampling and calculation of parasitism rate

Five to ten handfuls of about eight top branches of Cistus monspeliensis L. 1753 (cut at 50 cm below the end of the branch) were sampled in four localities (Fig. 1). These localities were part of a larger field survey of population dynamics of P. spumarius in Corsica. We targeted C. monspeliensis to maximize our chances to find eggs of P. spumarius. Indeed, we demonstrated in a previous study that, in Corsica, adults of P. spumarius seemed to be mainly associated with this species (Cruaud et al., 2018). Sampling was performed between the 12th and the 15th of February 2019. The back of each leaf (about 900 leaves per handful of branches) was inspected in the laboratory for whitish clusters, which were retained and inspected under a binocular microscope to confirm the presence of eggs of P. spumarius (Appendix S1; 109, 148, 167, 187 eggs obtained per site, hence a total of 611 eggs monitored). The morphological identification of P. spumarius eggs and first stage nymphs was performed using the descriptions of Weaver & King (1954) (Appendix S1). The pieces of leaf containing the eggs were placed on filter papers in unaerated Petri dishes (i.e., without spur) at room temperature (20.2 ± 1.5 °C), with natural light. Filter papers were kept moist by adding drops of water when necessary. Hatching was monitored every morning from the 18th of February to the 15th of March 2019. Emerging nymphs and parasitoids were killed and stored in 70% Ethanol at 4 °C. Parasitism rates were computed in each locality with the following formula: $Parasitismrate=\frac{Numberofparasitizedeggs}{Numberofparasitizedeggs+Numberofunparasitizedeggs}$ (Costello & Altieri, 1995).

Morphological identification of the parasitoids

Identification to species was performed using the Ooctonus keys by Triapitsyn (2010) and Huber (2012). Specimens were desiccated using HMDS (Heraty & Hawks, 1998) and glued on grey cards. Imaging was performed with a Keyence digital microscope (VHX-5000 Camera color CMOS and the VH-Z100UT lens). Images were then edited in Adobe Photoshop CS6©software.

Molecular identification of the parasitoids

Six individuals were used for molecular identification. Three of them were handled individually (sample codes = XMES00042_0101, XMES00077_0101, XMES00091_0101) and the remaining three were pooled to increase DNA yield (sample code = XMES00041_0189). Total genomic DNA was isolated using the Qiagen DNeasy Blood & Tissue kit without destruction of the specimens. We followed manufacturer’s protocol with the following modifications. Samples (whole insects, without dissection or crushing) were incubated overnight in an Eppendorf thermomixer (temperature = 56 °C, mixing frequency = 300 rpm). To increase DNA yield, two successive elutions (50 µL each) were performed with heated buffer AE (56 °C) and an incubation step of 15 min followed by centrifugation (6,000 g for 1 min at room temperature; see Cruaud et al. (2019) for a detailed description of the protocol). Eppendorf microtubes LoBind 1.5 ml were used for elution and to store DNA at minus 20 °C until PCR amplification. Vouchers were deposited at Centre de Biologie pour la Gestion des Populations (CBGP), Montferrier-sur-Lez, France. The mitochondrial Cytochrome c oxidase I standard barcode fragment (COI) was amplified with a cocktail of M13-tailed primers as detailed in Germain et al. (2013). Unpurified PCR products were sent to Eurofins MWG Operon (Ebersberg, Germany) for sequencing using the M13F and M13R primers (Germain et al., 2013; Ivanova et al., 2007). Both strands for each overlapping fragment were assembled in Geneious v11.1.4 (https://www.geneious.com). Geneious was also used to translate consensus sequences to amino acids to detect premature codon stops. All COI sequences available on BOLD (Ratnasingham & Hebert, 2007) for Ooctonus species were downloaded (last access July 12, 2019) and aligned with the newly generated sequences using MAFFT v7.245 (Katoh & Standley, 2013). A maximum likelihood tree was inferred with raxmlHPC-PTHREADS-AVX version 8.2.4 (Stamatakis, 2014). A rapid bootstrap search (100 replicates) followed by a thorough ML search (-m GTRGAMMA) was conducted. Tree visualization and annotation was performed with TreeGraph 2.13 (Stöver & Müller, 2010).

Species distribution modelling framework

Occurrences of O. vulgatus were retrieved from the literature and the GBIF database (GBIF.org, 2019) (Tables S1 and S2). Two hundred and five occurrences were obtained from the literature (Table S2), eight of which were not included in the analysis as no geographic coordinates were available. Forty occurrences were obtained from GBIF (last access: 22 August 2019; Table S2), but were all discarded due to dubious identification or lack of information on sample origin. Therefore, no occurrence obtained from GBIF could be included in the analysis.

We fitted a correlative model linking different climate descriptors to species occurrences. The Maxent algorithm was chosen to conduct analyses because it does not require absence data (i.e., locations in which we can presume that a species is truly absent) (Phillips, Anderson & Schapire, 2006). We summarized below the main step of our analysis and details are provided in Appendix S2. The mean temperature and precipitation of the wettest, driest, warmest, and coldest quarters as well as precipitation seasonality were extracted from the Worldclim 2.0 database (Fick & Hijmans, 2017) and used as bioclimatic descriptors (Hijmans et al., 2005). In absence of formal knowledge about climatic factors constraining O. vulgatus distribution, we constituted three sets of bioclimatic variables and performed modelling with each of them (Qiao, Soberón & Peterson, 2015; Godefroid et al., 2019). The first set (CLIM1) comprised the mean temperature of the wettest, driest, warmest, and coldest quarters to reflect the impact of temperature constraints on distribution. To highlight the precipitation constraint, we added the precipitation seasonality to CLIM1 and constituted the second set (CLIM2). Finally, we built a third set (CLIM3) by assembling CLIM1 and the precipitation of the wettest, driest, warmest, and coldest quarters to fully account for both extreme temperatures and precipitations in the species distribution models (SDMs). The Maxent algorithm requires a set of locations where the species has been found (here, a random 70% of the available occurrences, the other 30% being used for model validation) and a set of locations where no information about the presence of the species are available (referred to as background points). A total of 10,000 background points were randomly generated in North America and Europe. To render complex response to environmental constraints while reducing model overfitting we first fitted 48 Maxent models using six regularization multiplier (RM) combinations (L, LQ, H, LQH, LQHP, LQHPT with L = linear, Q = quadratic, H = hinge, P = product and T = threshold) and feature class (FC) values (eight values ranging from 0.5 to 4 with increments of 0.5) (Radosavljevic & Anderson, 2014). Optimal FC and RM combinations were determined for each of the three bioclimatic datasets (CLIM1–CLIM3) using the R language (R Core Team, 2019) and the package ENMeval (Muscarella et al., 2014). Optimal parameters were then used to fit a set of 10 replicate Maxent models using 70% of the dataset. The performance of each model was evaluated using the remaining 30% of occurrences using the area under the receiver–operator curve (AUC, Fielding & Bell, 1997) and the true skill statistics (TSS, Allouche, Tsoar & Kadmon, 2006). Models with AUC <0.8 were excluded from further analyses (Vicente et al., 2013). Habitat suitability maps (logistic output ranging from 0 to 1) were transformed into binary projections using the threshold that optimized the TSS statistics on the testing data (Guisan, Thuiller & Zimmermann, 2017). Maxent replicate models were fitted and evaluated using the R package biomod2 (Thuiller et al., 2009).

Two different outputs were generated using the set of model prediction. (i) Binary predictions were averaged to produce the committee (consensus) averaging (Araújo & New, 2007; Marmion et al., 2009) showing the likelihood of the presence of O. vulgatus. This consensus model ranges from 0 (all the models predict absence) to 100% (all the models predict presence) and (ii) the median of the logistic outputs (Guisan, Thuiller & Zimmermann, 2017) of the models that depicts the climate suitability across the different models.

Parasitism rates

Out of the 611 eggs monitored, 437 (i.e., 71.5%) hatched. 277 (63.4%) gave rise to P. spumarius nymphs and parasitoids emerged from 160 eggs (36.6%). All parasitoids were identified as O. vulgatus (Fig. 2). No parasitoid emerged from eggs collected in one of the four localities. We observed parasitism rates of 20.5, 48.9 and 69.0% in the three other localities (Fig. 1).

Guidelines for the identification of O. vulgatus

To help identification, we list below the main features that differentiate O. vulgatus from its closest relatives. The genus Ooctonus has been recently revised in the Palearctic and Nearctic regions respectively by Triapitsyn (2010) and Huber (2012). Ooctonus can be distinguished from other genera of Mymaridae by the following set of characters: tarsi 5-segmented, propodeum with diamond-shaped pattern of carinae (Fig. 2F), fore wing venation about one-third the wing length (Fig. 2C), with short marginal and stigmal vein, parastigma with hypochaeta next to proximal macrochaeta (Huber, 2012). In the Holarctic region, O. vulgatus can be distinguished from other species of Ooctonus by the following unique combination of features (Fig. 2): vertex without stemmaticum; mesoscutum without median groove; posterior part of scutellum and frenum smooth with weak sculpture laterally; metanotum and propodeum without reticulate sculpture; propodeum without median carina, but with a pentagonal areole formed by dorsolateral carinae; short petiole, 0.9–1.2x as long as metacoxa; forewing at least slightly truncate apically; females funicle with multiporous placoid sensilla (mps) on F7 and F8 only, F5 and F6 without mps; single row of six bullae inside the female clava; ovipositor at most 1. 4 × as long as metatibia and only slightly exerted beyond apex of gaster.

Molecular identification of the parasitoid

Barcode sequences were successfully generated from all samples. All sequences were identical. Phylogenetic analysis confirmed that the most likely identification was O. vulgatus (Fig. S1).

Species distribution modelling

A total of 200 occurrences (197 obtained from the literature plus the three localities where we sampled O. vulgatus) (Fig. 3A) were used to model the distribution of O. vulgatus in Europe. The optimal Maxent parameters were RM = 4 and FC = hinge; RM = 4 and FC = hinge and RM = 2.5 and FC = hinge for CLIM1, CLIM2, and CLIM3, respectively. With the exception of one model of CLIM2, all models based on these optimal values yielded AUC values >0.8, which indicated that the different bioclimatic data subsets performed well. The consensus model was therefore computed from a set of 29 estimates of climate suitability.

Figure 3B shows the median of the climate suitability values for the 29 models considered. Figure 3C depicts the proportion of the 29 models indicating that the climate is suitable for O. vulgatus. Both Figs. 3B and 3C show that the climate is favorable in very large areas covering most of Western Europe and around the Black Sea. These areas are overlapping with the geographical range of P. spumarius (Cruaud et al., 2018).