Preliminary Investigation of Species Diversity of Rice Hopper Parasitoids in Southeast Asia

Ongoing intensification of rice production systems in Southeast Asia is causing devastating yield losses each year due to rice hoppers. Their continuing development of immunity to resistant rice varieties and pesticide applications further complicates this problem. Hence, there is a high demand for biological control agents of rice hoppers. Egg parasitoid wasps are among the most important natural enemies of rice hoppers, such as Nilaparvata lugens and Nephotettix spp. However, our knowledge of their diversity is still very limited, due to their small size and the lack of available morphological information. Classifying these parasitoids is the first step to properly understanding their role in the rice agroecosystem. We used traditional morphological identification, as well as DNA sequencing of the 28S rRNA and the COI genes, to investigate the diversity of four important hopper egg parasitoid genera in the Philippines. Parasitoids of the genera Anagrus, Oligosita, Gonatocerus, and Paracentrobia were collected in eight study landscapes located in Luzon. Our findings illustrate that characterization of species diversity using morphological and molecular analyses were concordant only for the genus Paracentrobia. The genera Anagrus and Gonatocerus exhibited more genetic diversity than estimated with the morphological analysis, while the opposite was observed for Oligosita. This is the first study investigating the molecular diversity of rice hopper parasitoids in the Philippines. More research combining morphological, behavioral, and molecular methods, as well as the establishment of a comprehensive DNA database, are urgently needed to assess the performance and suitability of these organisms as biocontrol agents.


Introduction
Rice is the main food resource for more than half of the world's population [1,2]. This makes the rice production system of Asia one of the most important food production systems on Earth [3]. The brown planthopper (BPH; Nilaparvata lugens, Stål 1854) and the green leafhopper (GLH; Nephottetix spp.) are among the economically most important rice pests. These insects cause immense damage in the Asian rice paddies through xylem sap feeding, resulting in wilting, and subsequently, the death of the rice crops [4], as well as through the transfer of devastating viruses among fields [5,6].
So far, the introduction of rice varieties resistant to BPH and GLH has not been successful, as these pests rapidly adapt to the new varieties [7,8]. The misuse of pesticide applications further enhances the problem by disturbing the rice agroecosystem, which can increase the outbreak risk of hoppers [7,9,10].
As an alternative, integrated pest management (IPM) strategies can be used to maximize the pest's mortality by enhancing the role of their biological control agents. BPH and GLH are typically attacked by a wide range of natural enemies, such as spiders, predatory bugs, dragonflies, and egg parasitoid wasps from the families Mymaridae and Trichogrammatidae [11]. Among these enemies, parasitoid species are of particular interest. These organisms are very mobile, and can disperse over large distances [12]. The adults feed on pollen, nectar or the honeydew of their sap-feeding hosts, whereas their larvae develop in the hopper eggs and disrupt the hoppers' life cycles at the earliest possible stage [11]. In two previous studies on egg parasitoids of rice hoppers, egg parasitism levels of more than 60% were observed [13,14]. These findings are supported by Drechsler and Settele [15], who showed that parasitoids can play a major role in controlling hopper pests in rice agroecosystems.
Despite their importance, knowledge about species composition and diversity of parasitoid wasps in the tropics is still limited (see Nishida, Wongsiri and Wongsiri [16], as well as Gurr et al. [11]). Morphological species identification requires extensive taxonomic expertise, and is hampered by the small size (<1.5 mm) of the wasps as well as the limited amount of literature. To date, molecular information on rice hopper parasitoids from the Philippines is completely lacking. Molecular methods have become a promising tool to resolve species identities [17][18][19]. Different markers, such as the mitochondrial cytochrome c oxidase I (COI) and fragments of the small or large subunits of the ribosomal RNA genes (i.e., 18S or 28S rRNA genes), have been applied to construct hymenopteran and dipteran parasitoid phylogenies [20,21]. However, molecular analyses not only improve species identification, but may also reveal cryptic species [22,23]. Smith et al. [24] suggested combining barcoding with morphology and natural history. A similar conclusion was drawn by Padial and colleagues [25] in their review on an integrative taxonomy for the improvement of species discovery and description.
In a previous study, we investigated parasitoid wasps of Nilaparvata lugens (BPH) and Nephotettix spp. (GLH) in eight rice production landscapes located in Luzon, Philippines [26]. We found that BPH was parasitized by the Chalcid genera Anagrus and Oligosita, while GLH was parasitized by the chalcidoid genera Gonatocerus and Paracentrobia. In the present study, we analyzed the diversity of these genera with traditional morphological and molecular techniques. The present study provides the first molecular identification of parasitoid wasps in rice paddies in the Philippines, combined with a morphological identification based on dichotomous keys.

Study System and Sampling
This study was embedded in the project LEGATO, which focused on sustainable rice production [27]. Parasitoids were collected in eight study landscapes located in the Laguna province, Luzon, Philippines (Supplementary Figure S1). Sampling took place during the rice growing and fallow periods of the dry season from February to June 2013. In brief, rice plants of the variety Taichung Native (1) (TN1) were grown in a greenhouse for 6 weeks, trimmed to three tillers, and covered with small tubular insect cages (85 cm high, 15 cm diameter). Greenhouse cultures of BPH and GLH were reared on TN1, as previously described by Heinrichs et al. [28]. The BPH culture exclusively consisted of Nilaparvata lugens, while the GLH culture consisted of Nephotettix virescens (Distant) and Nephotettix nigropictus (Ståhl). Both hopper populations originated from wild individuals caught in the rice fields of the Laguna Province. Four gravid females of either BPH or GLH were released into each cage for 48 h to lay eggs. Subsequently, the plants were transferred to three different plots per study site to cover the naturally occurring landscape diversity in the Laguna province. Three plants infested with BPH eggs and three plants infested with GLH eggs were distributed randomly within each plot.
After 72 h, all plants were returned to the greenhouse and placed back inside separate insect cages. The parasitoids emerged 13-17 days after the field exposure of the eggs. Due to positive phototaxis behavior, parasitoids aggregated in small transparent glass vials, which were placed upside down on top of the sleeved cages. Adult parasitoids were collected daily within this period and immediately stored at −80 • C. The dead parasitoids were identified to genus under a stereo microscope using morphological keys [29][30][31]. In total, 19,455 parasitoids were collected and separated into the four genera Paracentrobia, Gonatocerus, Oligosita, and Anagrus [26].

Preliminary Morphological Identification
Fifty parasitoids from each of the four identified genera (Paracentrobia, Gonatocerus, Oligosita, and Anagrus) were randomly selected. The 200 parasitoids were permanently slide-mounted using Faure's medium, and identified to species level using a microscope. A specialist at the International Rice Research Institute (IRRI) helped with the identification using the key published in Heong and Hardy [30].

DNA Extraction and Amplification
Total genomic DNA was extracted from 267 newly selected whole single female parasitoid individuals (Supplementary Table S2) employing the Phire Animal Tissue Direct PCR Kit (Thermo Scientific, Waltham, MA, USA). Parasitoids were suspended in 20 µL TE buffer (100 mM Tris, 10 mM EDTA). To increase DNA extraction efficiency, 1 µL proteinase K (20 mg mL −1 ) was added. Parasitoid samples were carefully homogenized with sterile micro pestles, and incubated at room temperature overnight.
Two independent gene fragments were amplified from the extracted DNA: one located on the COI subunit I (COI I, amplicon length around 670 bp) and the other one on the expansion regions D2-D3 of the 28S ribosomal subunit (28S-D2/28S-D3, amplicon length approximately 610 bp). The COI region was amplified using the primer pair HCO2198/LCO1490 [32]. The PCR reaction mixture (25 µL) contained 2.5 µL of 10-fold Ex Taq Buffer (Takara Biotechnology, co., LTD, Dalian, China), 25 mM MgCl 2 , 2.5 mM of each of the four dNTPs (deoxynucleotide triphosphates), 10 µM of each primer, 1 U TaKaRa Ex Taq polymerase (Takara Biotechnology), and approximately 25 ng of parasitoid DNA. The following thermal cycling scheme was used: initial denaturation at 94 • C for 3 min, 5 cycles of denaturation at 94 • C for 45 s, annealing at 45 • C for 45 s, followed by elongation at 72 • C for 1 min, and 25 cycles of denaturation at 94 • C for 45 s, annealing at 50 • C for 45 s, followed by another elongation at 72 • C for 1 min. The final extension was carried out at 72 • C for 5 min.
The D2-D3 region of the 28S rRNA gene was amplified as described for the PCR above, using the primers D2-3549 [33] and D2-4068 [34]. The following thermal cycling scheme was used: initial denaturation at 94 • C for 3 min, 30 cycles of denaturation at 94 • C for 45 s, annealing at 58 • C for 45 s, followed by elongation at 72 • C for 1 min. The final extension was carried out at 72 • C for 6 min. Negative controls were performed by using the reaction mixture without template. PCR products were controlled for appropriate size and then purified using the peqGOLD Gel Extraction kit as recommended by the manufacturer (Peqlab, Erlangen, Germany; now VWR).
Sequencing was performed at the Göttingen Genomics Laboratory using an ABI 3730xl system and BigDye terminator chemistry version 3.1 (Thermo Fisher Scientific).

Data Analysis
Quality control and trimming of the forward and reverse DNA sequences were performed with Gap5 v 1.2.14-r [35]. DNA sequence alignment was performed using MEGA 6 [36]. Additional Chalcidoidea (the super family containing all four genera) sequences were used as reference material. They were retrieved by using the National Centre for Biotechnology Information (NCBI) BLAST tool [37] and the BOLD Identification System (http://www.boldsystems.org) (Supplementary Table S3). The best model to calculate parasitoid phylogeny was determined assuming partial deletion, site coverage cut-off of 95%, and the branch swap filter set to "very strong". Maximum likelihood trees were generated under the assumption of the best fitting model with 1000 bootstraps. Finally, median-joining networks were constructed using the software NETWORK v. 4.6.1.2 [38].
The genetic diversity between and within species was estimated using the reference sequences from the superfamily Chalcidoidea (Supplementary Table S3). The appropriate models for the calculation of the pairwise molecular distances between individuals were determined with MEGA 6. The 28S rRNA distances were calculated based on the K2P + G model. COI distances were calculated based on the T92 + G model. The values were compared using a Kruskal-Wallis test, followed by a Dunn's post hoc test [39] in R version 3.2.3 [40].
To estimate, whether the genetic diversity of the parasitoids was covered by the sample size used in the analyses, we performed a rarefaction analysis (Supplementary Figures S4 and S5) in R version 3.2.3 using the packages vegan [41] and drc [42].

Molecular Analysis
Of the 267 specimens used for the molecular analysis, we successfully sequenced a total of 105 parasitoid samples (Supplementary Table S2). We failed to sequence 162 specimens, due to the poor quality of the extracted DNA, or a failure to extract DNA from the samples. A total of 86 (COI) and 105 (28S) sequences were used to create the maximum likelihood trees. The final dataset for comparison of the two gene fragments included the sequences from 74 parasitoid individuals.
Rarefaction analysis for the local hopper egg parasitoids genetic diversity revealed that the majority of the 28S rRNA gene haplotypes was recovered by the surveying effort (Supplementary Figure S2). By contrast, the rarefaction curve for the COI gene haplotypes was not saturated (Supplementary Figure S3).

Comparison between the Morphological and Molecular Approach
The molecular analyses based on COI and 28S rRNA genes revealed that the sequences clearly segregated according to the morphologically pre-assigned genera (Figures 1-4, Supplementary  Figures S4 and S5). Oligosita exhibited the lowest genetic diversity, with three highly similar sequences in the COI gene analysis and only one haplotype in the 28S rRNA gene analysis. In contrast to that, three Oligosita species were identified using the morphological approach.
Paracentrobia sequences varied only slightly, with two highly similar haplotypes in the 28S rRNA gene analysis and three highly similar haplotypes in the COI gene analysis (Figures 1-4, Supplementary Figures S4 and S5). This result indicates that all samples of Paracentrobia that we analyzed morphologically and molecularly belonged to a single species.
The sequences for both genes from Gonatocerus exhibited more variability compared to the morphological data alone. Three clades were unambiguously identified, with one cluster occurring prevalently (88.2% for the 28S gene, 83.3% for the COI gene) compared to the other two clusters (Figures 1-4, Supplementary Figures S4 and S5). Anagrus exhibited more variability, and was the most genetically diverse genus: 7 and 12 different haplotypes were identified in the 28S rRNA gene analysis and in the COI gene analysis, respectively (Figures 1 and 2 sequences in the COI gene analysis and only one haplotype in the 28S rRNA gene analysis. In contrast to that, three Oligosita species were identified using the morphological approach. Paracentrobia sequences varied only slightly, with two highly similar haplotypes in the 28S rRNA gene analysis and three highly similar haplotypes in the COI gene analysis (Figures 1-4, Supplementary Figures S4 and S5). This result indicates that all samples of Paracentrobia that we analyzed morphologically and molecularly belonged to a single species.
The sequences for both genes from Gonatocerus exhibited more variability compared to the morphological data alone. Three clades were unambiguously identified, with one cluster occurring prevalently (88.2% for the 28S gene, 83.3% for the COI gene) compared to the other two clusters (Figures 1-4, Supplementary Figures S4 and S5). Anagrus exhibited more variability, and was the most genetically diverse genus: 7 and 12 different haplotypes were identified in the 28S rRNA gene analysis and in the COI gene analysis, respectively (Figures 1 and 2          The degrees of genetic diversity found in the genera Anagrus and Gonatocerus sequences were particularly high. Within group pairwise distances were 0.075 ± 0.012 SE (28S rRNA gene) and 0.050 ± 0.009 SE (COI gene) for Anagrus, and 0.069 ± 0.012 SE (28S rRNA gene) and 0.100 ± 0.016 SE (COI gene) for Gonatocerus. By contrast, within group pairwise distances were 0.000 ± 0.000 (28S rRNA gene) and 0.005 ± 0.003 SE (COI gene) for Oligosita, and 0.001 ± 0.001 SE (28S rRNA gene) and 0.009 ± 0.004 SE (COI gene) for Paracentrobia ( Figure 5). Paracentrobia and Oligosita were significantly different from the congeneric Chalcidoidea, but not from the conspecific Chalcidoidea ( Figure 5, Table 2). The genetic differences within Gonatocerus and Anagrus were in the same order of magnitude as the congeneric Chalcidoidea data.  The degrees of genetic diversity found in the genera Anagrus and Gonatocerus sequences were particularly high. Within group pairwise distances were 0.075 ± 0.012 SE (28S rRNA gene) and 0.050 ± 0.009 SE (COI gene) for Anagrus, and 0.069 ± 0.012 SE (28S rRNA gene) and 0.100 ± 0.016 SE (COI gene) for Gonatocerus. By contrast, within group pairwise distances were 0.000 ± 0.000 (28S rRNA gene) and 0.005 ± 0.003 SE (COI gene) for Oligosita, and 0.001 ± 0.001 SE (28S rRNA gene) and 0.009 ± 0.004 SE (COI gene) for Paracentrobia ( Figure 5). Paracentrobia and Oligosita were significantly different from the congeneric Chalcidoidea, but not from the conspecific Chalcidoidea ( Figure 5, Table 2). The genetic differences within Gonatocerus and Anagrus were in the same order of magnitude as the congeneric Chalcidoidea data. The degrees of genetic diversity found in the genera Anagrus and Gonatocerus sequences were particularly high. Within group pairwise distances were 0.075 ± 0.012 SE (28S rRNA gene) and 0.050 ± 0.009 SE (COI gene) for Anagrus, and 0.069 ± 0.012 SE (28S rRNA gene) and 0.100 ± 0.016 SE (COI gene) for Gonatocerus. By contrast, within group pairwise distances were 0.000 ± 0.000 (28S rRNA gene) and 0.005 ± 0.003 SE (COI gene) for Oligosita, and 0.001 ± 0.001 SE (28S rRNA gene) and 0.009 ± 0.004 SE (COI gene) for Paracentrobia ( Figure 5). Paracentrobia and Oligosita were significantly different from the congeneric Chalcidoidea, but not from the conspecific Chalcidoidea ( Figure 5, Table 2). The genetic differences within Gonatocerus and Anagrus were in the same order of magnitude as the congeneric Chalcidoidea data.

Discussion
The molecular analyses of the diversity within the four Chalcidoidea genera studied showed discrepancies with the morphological analysis. This finding is in line with previous studies on hymenopteran parasitoids [23,43,44]. For example, Mottern and Heraty [45] found that one species of the parasitoid previously described as Cales noacki were actually ten different Cales species. Although previous studies also showed that DNA barcoding can be misleading, resulting in potentially more species than actually exist [46,47], the COI gene, as the classical barcoding gene, has shown its utility in many studies [43,48,49], and constitutes the basis for the BOLD database (www.barcodinglife.org). The 28S rRNA gene is valuable to distinguish among closely related species of Chalcidoidea [50][51][52]. In addition, the 28S rRNA gene is generally more conserved than the COI gene for this group of insects [44,45,53]. This observation is supported by our results, as we found that sampling saturation was only reached for the 28S rRNA gene but not for the COI gene.
Once a reliable and comprehensive DNA database is established, the identification of hopper parasitoids using DNA sequencing will be a vital step towards assessing the performance and suitability of these organisms as biocontrol agents [54,55]. Our study shows that an accurate identification relies on both molecular and morphological techniques. Similar conclusions were drawn previously [21,23,56].
For instance, to date, Anagrus species are described as generalist parasitoids that can switch between alternative hopper species in rice agroecosystems [57,58]. Yet, the strong diversity found in the present study by using molecular tools suggests that the Anagrus genus is a complex of cryptic species that may be host specialists. This distinction is of high importance to control the hoppers effectively. In addition, the discovery of unknown cryptic species could expand the list of potential biological control agents [45].
Although we were not able to validate the existence of new species of Anagrus and Gonatocerus by mating tests and/or further morphometric analyses, our data strongly indicates that the two genera include unknown species. The level of within-genus genetic difference, for both 28S rRNA and COI gene sequences among samples from the Anagrus and Gonatocerus genera, exceeded the threshold of 0.17-2% within-sequence variation, which is generally accepted to delineate individuals of the same species [17,49,59]. According to the literature [17,49,59], genera from the Chalcidoidea superfamily diverge by 5.8-11.25%. Thus, the level of genetic differences we observed for Anagrus and Gonatocerus suggests that these two genera were composed of individuals from different species. When comparing the genetic analyses to the morphological results, we found strong evidence that there was at least one species of Anagrus and Gonatocerus that was not detected by the morphological analysis. Interestingly, we observed the opposite for Oligosita samples. In the latter genus, it is possible that one or more species were not detected due to the low success rate of DNA amplification. It is also possible that we observed a discrepancy between morphological and molecular analyses, because they were performed on a different set of specimens. As one Oligosita species was much more prevalent than the other two in the morphological dataset, we cannot exclude that we analyzed the most common species only, in the molecular analysis.

Conclusions
In conclusion, our results clearly demonstrate that molecular identification should be used in combination with morphological methods and mating tests for assessing the diversity of rice hopper parasitoids. However, further studies using an integrative approach are needed to cover the whole diversity of parasitoids, as well as to find sustainable solutions to problems caused by BPH and GLH. To validate the potential application of parasitoid wasps as biocontrol agents, it is of crucial importance to have a comprehensive knowledge on their ecology and diversity. We hope that this study will encourage further research by providing the first barcodes for egg parasitoid species from rice paddies in Southeast Asia.
Supplementary Materials: The following are available online at http://www.mdpi.com/2075-4450/9/1/19/s1, Figure S1: Distribution of the eight study sites throughout the Laguna province, Figure S2: Rarefaction curve for the 28S rRNA gene, Figure S3: Rarefaction curve for the COI gene, Figure S4: Maximum likelihood tree for the 28S rRNA sequences, Figure S5: Maximum likelihood tree for the COI sequences, Table S1: Morphometric results for parasitoids throughout the Laguna province, Philippines, Table S2: Sequences obtained from parasitoids of the rice fields in the Philippines, Table S3: Sequences retrieved from the NCBI GenBank.